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ABSTRACT 

We present the first observations of H 13 CN(1 - 0), H 13 CO+(l - 0) and SiO(2 - 1) in NGC6240, 
obtained with the IRAM PdBI. Combining a Markov Chain Monte Carlo (MCMC) code with Large 
Velocity Gradient (LVG) modelling, and with additional data from the literature, we simultaneously 
fit three gas phases and six molecular species to constrain the physical condition of the molecular 
gas, including mass—luminosity conversion factors. We find ~ 10 1o Mq of dense molecular gas in 
cold, dense clouds (Ik ~ 10 K, uh 2 ~ 10 6 cm -3 ) with a volume filling factor < 0.002, embedded in a 
shock heated molecular medium (Tk ~ 2000K, tih 2 ~ 10 3 6 cm -3 ), both surrounded by an extended 
diffuse phase (T k ~ 200K, n H2 ~ 10 2 5 cm _3 ). We derive a global aco = 1.5,;} with gas masses 
log 10 (M/[Af 0 ]) = IO.I^q q, dominated by the dense gas. We also find «hcn = 32f§, which traces 
the cold, dense gas. The [ 12 C]/[ 13 C] ratio is only slightly elevated (98§ 30 ), contrary to the very high 
[CO]/[ 13 CO] ratio (300-500) reported in the literature. However, we find very high [HCN]/[H 13 CN] 
and [HCO + ]/[H 13 CO + ] abundance ratios (3002 qq) which we attribute to isotope fractionation in the 
cold, dense clouds. 


1. INTRODUCTION 

The term “dense gas” typically refers to molecular gas 
with nn 2 > 10 4 cm -3 , which is the typical density 
of molecular cloud cores bound by their self-gravity 
(e.g., Lada 1992). This phase offers significant shielding 
against photo-dissociation while being heated and par¬ 
tially ionised by cosmic rays, leading to chemistry dis¬ 
tinct from elsewhere in the ISM. It is best traced by the 
rotational transitions of molecules such as HCN, HCO + , 
and CS, which can reach high abundances relative to H 2 
(10 -8 —10~ 6 ) under these conditions (Nguyen et al. 1992; 
Solomon et al. 1992; Boger & Sternberg 2005). Due to 
both high densities and line trapping the low—J transi¬ 
tions of these lines (in particular HCN) are almost ex¬ 
clusively optically thick, so that approaches to estimat¬ 
ing gas mass developed for CO (Young & Scoville 1982; 
Papadopoulos et al. 2012; Bolatto et al. 2013) become 
applicable, with appropriate calibrations, allowing us to 
trace specifically dense molecular gas. 

Single-dish surveys of primarily the HCN, but also the 
HCO + , J = 1 — 0 line towards samples of local (ul- 
tra)luminous infrared galaxies ((U)LIRGs) suggest that 
a larger fraction of the molecular gas in these galaxies is 
in a dense state compared to that of normal star form¬ 
ing galaxies (Solomon et al. 1992; Gao & Solomon 2004; 
Garcfa-Burillo et al. 2012), where the HCN/CO J = 1 — 0 
luminosity ratio is adopted as a proxy of the dense gas 
mass fraction (typically ~ 0.1 in LIRGs/ULIRGs cf. 
~ 0.03 in normal star forming galaxies, Gao & Solomon 
2004; Greve et al. 2014). However, to date, only a liand- 
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ful of sources have been analysed using multi-line data¬ 
sets of dense gas tracers, in conjunction with radiative 
transfer analyses, to constrain the dense gas mass (Krips 
et al. 2008; Greve et al. 2009; Papadopoulos et al. 2014). 

Gas mass conversion factors typically assume: 1) 
something about the kinematical state of the gas; usu¬ 
ally that it resides within virialized clouds (Bolatto et al. 
2013, although see Papadopoulos et al. 2012 for a for¬ 
mulation that allows for a generalised kinematic state): 
perhaps not unreasonable for dense gas tracers but cer¬ 
tainly questionably for CO. 2) that the lines are opti¬ 
cally thick, which is generally a good assumption, at least 
for the J = 1 — 0 lines. 3) metallicity and/or molecu¬ 
lar abundances. These are particularly uncertain, be¬ 
ing strongly affected by the environment the emission 
is tracing, e.g., photon dominated regions (PDRs), X- 
ray dominated regions (XDRs), hot-cores or cold, dark 
clouds, and are usually a subject of interest in them¬ 
selves. There are also possible contaminations from other 
sources, such as mid-IR pumping of HCN (and perhaps 
HNC and HCO + ), which increases the low—J HCN line 
intensity (Aalto et al. 2007, 2012, 2015). As we move to 
higher redshifts the effect of the CMB becomes notice¬ 
ably non-linear, with the greatest effects on the low—J 
CO lines, and even in the local universe bright contin¬ 
uum backgrounds could be having significant effects on 
observed lines and line ratios (Papadopoulos et al. 2000, 
2010). HCO + is subject to a complicated network of ion 
chemistry and photochemical reactions, making it par¬ 
ticularly complicated molecular tracer (Viti et al. 2002; 
Papadopoulos 2007). 

Optically thin dense gas tracers, in particular isotopo- 
logues of the main species (e.g., H 13 CN and H 13 CO + ), 
can be especially powerful as they provide much tighter 
constraints on large velocity gradient (LVG) models 
when combined with their 12 C isotopologues. Although 
usually 40 — 100 x, and perhaps up to ~ 500 x, less abun¬ 
dant than their 12 C isotopologues, optical depth effects 
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mean that these lines are usually only ~ 10—15 x fainter, 
and are detectable with the IR.AM Plateau de Bure Inter¬ 
ferometer (PdBI) and certainly with the Atacama Large 
Millimeter/submillimeter Array (ALMA) for the closest 
LIRGs and ULIRGs. 

1.1. Chemical Modelling and Isotopologue Ratios 

If it can be accurately measured, the elemental 
[ 12 C]/[ 13 C] abundance ratio can be a powerful diagnos¬ 
tic of the evolutionary state of a galaxy, with a high 
[ 12 C]/[ 13 C] being evidence for a young, unprocessed ISM 
(Milam et al. 2005; Henkel et al. 2010). Recently, it 
has been suggested by Papadopoulos et al. (2011) that 
a high [ 12 C]/[ 13 C] ratio might be evidence for a cosmic 
ray dominated star formation paradigm, leading to a top 
heavy stellar initial mass function and thereby an ex¬ 
cess of massive, 12 C producing, stars. Typical values of 
[ 12 C]/[ 13 C] in the local universe are ~ 25 in the cen¬ 
tral Milky Way, 68 in the Milky Way as a whole and up 
to ~ 100 in local starbursts and > 100 in high redshift 
ULIRGs (Milam et al. 2005; Henkel et al. 2010, 2014). 
Typically, these ratios are measured using a combination 
of observations of CN, 13 CN, CO, 13 CO and C 18 0 lines. 
However, this measurement is non-trivial due to chemical 
effects leading to isotope fractionation under a variety of 
circumstances, as well as optical depth effects. 

The theory of 12 C/ 13 C isotope fractionation rests heav¬ 
ily, although not exclusively, upon isotope charge ex¬ 
change (ICE) reactions, such as: 

13 C+ + CO ^ C+ + 13 CO + 35 K, (1) 

13 CO + HCO+^ CO + H 13 CO++ 17K, (2) 

13 C + + CN ^ C + + 13 CN + 31K, (3) 

(Watson et al. 1976; Langer et al. 1978, 1984; Roueff 
et al. 2015). Since these forward reactions are slightly 
exothermic, in cold molecular gas these reactions can lead 
isotope fractionation. However, there are large and com¬ 
plex chemical networks associated with and linking these 
and other molecules, necessitating sophisticated chemi¬ 
cal modelling. There is also competition with selective 
photo-dissociation (SPD): for CO UV photo-dissociation 
is a line process (van Dishoeck & Black 1988), allow¬ 
ing for isotopologue specific self-shielding and thereby 
fractionating CO in the opposite direction to ICE. In 
contrast, CN photo-dissociation is a continuum process 
(Lavendy et al. 1987). 

Langer et al. (1984) presented time dependent chem¬ 
ical models of molecular cloud chemistry for a range of 
densities (5xl0 2 — 1x10 s cm -3 ), temperatures (6—80K), 
metallicities and C/O abundances, with the aim of mod¬ 
elling 13 C fractionation. They found large variations 
between species and with conditions, as well as multi¬ 
ple degeneracies. Three distinctly behaved carbon pools 
emerged: CO, HCO + , and other carbon bearing species 
such as C + , H 2 CO, HCN, CS, CH etc., with each pool 
broadly following a given trend in fractionation. 

Their key findings were the confirmation that a small 
decrease in the [CO]/[ 13 CO] ratio from ICE can, due 
to the high relative abundance of CO, rapidly deplete 
the pool of available 13 C, leading to very elevated ratios 
such as [HCN]/[H 13 CN], which attained ratios as high 


as 6x that of [CO]/[ 13 CO]. HCO+ on the other hand 
exists in a complex equilibrium with both other pools 
and can be fractionated in either direction, although to 
a lesser extent. It is also very susceptible to variations 
in metallicity and the C/O ratio, but was never found 
to deviate from the [CO]/[ 13 CO] ratio by more than a 
factor of 1.5. Combining their results, they predicted a 
bracketing of the true [ 12 C]/[ 13 C] ratio: 

[CO]/[ 13 CO] < [ 12 C]/[ 13 C] < [HCN]/[H 13 CN], (4) 

while [HCO + ]/[H 13 CO + ] is impossible to interpret in iso¬ 
lation, and extremely challenging even when combined 
with observations of multiple species and their isotopo¬ 
logue ratios. 

Observational studies over the subsequent three decades 
have presented an even more complex picture. Milam 
et al. (2005) observed CN and 13 CN along lines of sight 
to dense clouds in the Milky Way, finding no difference 
in the isotopologue ratios of CN, CO and H 2 CO, instead 
finding that all three ratios increase from ~ 25 in the 
Galactic Centre to ~ 130 at a Galactic radius of 16.4 kpc, 
ostensibly due to the ISM being more evolved (and hence 
enriched in 13 C) near the Galactic centre. 

On the other hand, Sheffer et al. (2007) used HST 
UV observations to measure [CO]/[ 13 CO] over 25 lines 
of sight, finding evidence of CO isotope fractionation in 
both positive and negative directions, as well as find¬ 
ing that the [CN]/[ 13 CN] ratio appeared to anti-correlate 
with [CO]/[ 13 CO], as predicted by Langer et al. (1984). 

Similarly, Ritchey et al. (2011) found significant ev¬ 
idence for anti-correlated fractionation of CO and CN 
along lines of sight to diffuse molecular clouds, such as 
£ Oph, which presents [CO]/[ 13 CO] = 167 ± 15 and 
[CN]/[ 13 CN] = 47.3+44 (Lambert et al. 1994; Crane 
& Hcgyi 1988). Ritchey et al. (2011) identified CH+ 
as a potential candidate to reliably constrain the true 
[ 12 C]/[ 13 C] ratio as it is produced only in non-thermal 
processes 5 , so is insusceptible to fractionation. Indeed, 
towards £ Oph [CH + ]/[ 13 CH + ] = 67.5±4.5, c.f. the ISM 
standard of 68 (Crane et al. 1991; Milam et al. 2005). 

The most recent, state-of-the-art chemical modelling 
study of isotope fractionation was conducted by Roueff 
et al. (2015), who aimed to refine the understanding 
of [ 14 N]/[ 15 N] fractionation by cementing the related 
[ 12 C]/[ 13 C] fractionation chemistry. Their detailed chem¬ 
ical model included some of the most up-to-date reac¬ 
tions and chemical networks and explored time depen¬ 
dent isotope fractionation in cold (10 K), dense (2 x 10 4 
and 2 x 10 5 cm -3 ) molecular clouds, with an elemental 
[ 12 C]/[ 13 C] = 68. One of the significant changes is the 
addition of the reaction: 

HNC + C -> HCN ± C, (6) 

5 CH+ is formed via the reaction: 

C+ + H 2 -> CH+ + H, (5) 

which is endothermic with AE/fcg = 4640 K, so requires localised 
magnetohydrodynamic shocks to form in molecular clouds (Elitzur 
& Watson 1978, 1980; Ritchey et al. 2011) 
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(Loison et al. 2014), which couples HCN back into the 
chemical network post-production, allowing it to be frac¬ 
tionated at any time that the 13 C pool is depleted. In 
addition, Roueff et al. (2015) use an updated AE /ks for 
Equation 2 of 17 K, c.f. 9K in Langer et al. (1984) (Lohr 
1998; Mladenovic & Roueff 2014). 

The models of Roueff et al. (2015) have particularly 
interesting results for the isotopologue ratios of CN, 
HCN and HCO + . Like Langer et al. (1984), Roueff 
et al. (2015) found that in steady state [CO]/[ 13 CO] 
is only very slightly less than [ 12 C]/[ 13 C], at 67.4 and 
68 in the lower and higher density models respectively. 
CN fractionates up in the first model, and down in 
the second. HCN on the other hand is elevated in all 
of their models, with steady state ratios consistently 
2 ~ 3x [CO]/[ 13 CO], HCO + tends to equilibrate with 
[CO]/[ 13 CO] in all of the models. They also found that 
the H 2 ortho-para ratio has very significant effects on 
some species, such as CN, but not others, such as HCN. 

The fundamental conclusion of the observational and the¬ 
oretical studies is that isotope fractionation is possible, 
but that it is complex and hard to constrain due to the 
multitude of variables driving and opposing fractiona¬ 
tion. Even state-of-the-art chemical models are limited 
in usefulness unless they have been run for the specific 
cloud conditions in question due to the strong sensitiv¬ 
ity to temperature, density and metallicity. However, 
HCN appears to be reliably and strongly fractionated 
away from the elemental abundances, while HCO + can 
be fractionated in either direction and is very sensitive 
to the evolutionary stage of the molecular cloud. 

It is standard practice, when observing H 13 CN, to as¬ 
sume a canonical [HCN]/[H 13 CN] abundance ratio (e.g., 
Pirogov et al. 1995; Nakajima et al. 2011; Tunnard et al. 
2015). In light of the chemical modelling discussed above 
this approach appears to risk being significantly flawed: 
since HCN preferentially traces colder, denser gas, it is 
especially sensitive to the effects of isotope fractionation. 
Also, there is a degeneracy between the kinetic temper¬ 
ature and the isotopologue ratio for HCN/H 13 CN line 
ratios; assuming an incorrect [HCN]/[H 13 CN] runs the 
risk of significantly biasing derived HCN column densi¬ 
ties and estimates of the local kinetic temperature. 

The chemical studies also suggest that any claims of 
elevated elemental [ 12 C]/[ 13 C] ratios, especially in extra- 
galactic sources, must be interpreted with caution, and 
within not only the appropriate physical context, but 
also the appropriate chemical context (e.g., as in Henkel 
et al. 2014). There is also the potential for strong optical 
depth effects (Szucs et al. 2014), although that particu¬ 
lar challenge appears intractable with unresolved, single 
dish extragalactic observations. 

1.2. NGC6240 

NGC 6240 is an early stage major merger galaxy in the 
constellation of Ophiuchus. Hosting two active super- 
massive black holes and possessing an infra-red lumi¬ 
nosity As—looo/xm = 10 11 ' 73 L 0 it is a luminous infrared 
galaxy (LIRG) and a LINER galaxy (low-ionization nu¬ 
clear emission-line region) (Komossa et al. 2003; Lutz 
et al. 2003; Sanders et al. 2003; Veron-Cetty & Veron 


2006). Studies by Tacconi et al. (1999) and Tecza et al. 
(2000) have shown that while star formation is con¬ 
centrated in two nuclei, the molecular gas is predomi¬ 
nantly in a turbulent, thick inter-nuclear disk. They also 
demonstrated that the extremely luminous infrared H 2 
5(0) and 5(1) emission (L(H 2 ) ~ 2 x 1O 9 L 0 ) is most 
likely due to slow C shocks between gas clouds in this 
dense, turbulent disk; a finding confirmed by the study 
of the CO spectral line energy distribution (SLED), com¬ 
plete up to and including J = 13— 12, by Meijerink et al. 
(2013). 

As one of the closest LIRGs the molecular emission 
in NGC 6240 has been extensively studied. The find¬ 
ing of Greve et al. (2009) that the CO SLED required at 
least a two component fit has been extended upon by Pa- 
padopoulos et al. (2014), who fitted the CO SLED of CO 
from J = 1 — 0 to J = 13—12 with an “inside out” decom¬ 
position, where the results of HCN and HCO + J = 1 — 0 
to J = 4 — 3 SLED fits are used to constrain the high—J 
CO lines, and a third phase fitted to the remaining CO 
lines. Two of their key findings for NGC 6240 were a very 
high [ 12 C]/[ 13 C] ratio of 300 — 500 and that most of the 
H 2 gas mass is inconsistent with self-gravitating and pho- 
toelectrically heated gas, requiring either cosmic ray or 
turbulent heating. Indeed, Meijerink et al. (2013) model 
the CO SLED as being entirely due to shock heated gas. 
An independent two phase analysis was also conducted 
by Kamenetzky et al. (2014), who used a nested sam¬ 
pling routine to model only the CO SLED, finding a hot, 
diffuse phase and a cold, dense phase. 

Papadopoulos et al. (2014) argued that the high 
[ 12 C]/[ 13 C] (300-500) abundance ratio they found im¬ 
plies the existence of a very unusual environment in 
NGC 6240. They suggest that it is either a sign of a 
prolonged period of star formation with a top heavy ini¬ 
tial mass function (IMF) or an inflow of unprocessed gas 
from the outer edges of the galaxy. However, the massive, 
wide outflow in CO observed by Feruglio et al. (2013) in 
NGC 6240 suggests that an extensive inflow of gas would 
be difficult to maintain. That being said, Garcfa-Burillo 
et al. (2014) observed simultaneous inflows and outflows 
into/from the circunmuclear disk of NGC 1068, so this 
possibility cannot be entirely ruled out. 

As well as presenting a striking CO SLED, NGC 6240 
was shown to possess massive (~ 120 M 0 yr _1 ) molecular 
outflows at —600, +800kms _1 by Feruglio et al. (2013) 
in CO(l — 0), while Wang et al. (2014) found that the 
extended (5kpc) hard X-ray spectrum (~ 6keV) is due 
to ~ 2200kms _1 outflows, which they ascribe to the 
nuclear starburst with a supernova rate of ~ 2 yr _1 . 

This combination of extreme conditions make 
NGC 6240 a unique laboratory for star formation 
physics under unusual conditions. 

This paper is outlined as follows: we describe our obser¬ 
vations and their analysis in Section 2 and the line extrac¬ 
tion in Section 3. We introduce our Monte Carlo Markov 
Chain - Large Velocity Gradient (MCMC LVG) model 
and results in Section 4. The results are discussed in Sec¬ 
tion 5 and we outline our conclusions in Section 6. We 
include further details of our model in the Appendix. We 
adopt cosmological parameters h = 0.70, D m = 0.3 and 
Da = 0.7, giving a luminosity distance of 107 Mpc and 
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an angular scale of 0.494 kpc/" at 2 = 0.0245 (Wright 
2006; Papadopoulos et al. 2014). To avoid confusion, we 
note that throughout this work, unless explicitly stated 
otherwise, we present results as x%, where a and b de¬ 
note the values at the upper and lower limits on the 68% 
credible interval. 

2. OBSERVATIONS 

We obtained observations of H 13 CN(1 — 0), H 13 CO + (l — 
0) and SiO(2 — 1) in NGC6240 with the Plateau de 
Bure Interferometer (PdBI) on 2008 February 21 and 
2008 March 7 in the most extended configuration (PI: 
J. Gracia-Carpio), obtaining minimum and maximum 
baselines of 59.7 m (16.8 kA) and 760 m (215.8 kA) respec¬ 
tively, centred on the sky frequency of 84.6 GHz. Using 
natural weighting to maximise sensitivity we obtained a 
synthesised beam 1.86" x 0.74". 

The data were reduced in CLIC, part of the IRAM 
software package GILD AS 6 7 , while for the UV analysis 
and imaging we used the GILDAS program, MAPPING. 

3. SPECTRUM EXTRACTION AND DATA ANALYSIS 

Due to the poor signal to noise of the UV visibilities we 
analysed the data independently in both the UV plane 
and in the final imaged datacube. These two methods 
are outlined below, and produced spectra consistent on 
a channel by channel basis, within the channel uncertain¬ 
ties. 

Elliptical Gaussian and disk models were fitted to the 
UV visibilities on a channel-by-channel basis, using the 
GILDAS task UVJTT and 30 MHz channels. The fitted 
fluxes and fit rms’s were identical for the two models 
and we adopted the elliptical Gaussian model. Fitting 
of channels narrower than 30 MHz was attempted but 
produced pathological fits in multiple channels due to 
the poor S/N of the visibilities. The spectrum obtained 
is shown in Figure 1, where we fit two Gaussian profiles 
and a 0 th order polynomial baseline'. The SiO(2 — 1) 
line is clearly detected and there is a possible detection 
of H 13 CN(1 — 0). The significance of this detection is 
quantified in Section 3.2 below. Channel errors are shown 
at the Tier level and are unique for each channel, based 
on uncertainty in the fitted flux. The median channel 
error is ±0.8 mJy. 

A concern with the UV fit to each channel is 
the possibility of fitting to a different spatial centre, 
and thereby unintentionally applying a spatial pseudo¬ 
averaging across the channels. The majority of the chan¬ 
nels fitted a centre consistent within the uncertainties 
(~ ±0.06"), while the full range is 0.3" in Dec and 0.1" 
in RA, dominated by four outlying channels. These off¬ 
sets are < 25% of the fitted FWHMs, and < 5% of the 
aperture used for extracting the spectrum from the im¬ 
aged spectrum (see below), and so this is not a significant 
effect here. 

6 http://www.iram.fr/IRAMFR/GILDAS/. 

7 All reasonable fits with a 1 st order polynomial baseline had a 
gradient of almost zero. However, due to the randomised starting 
conditions we use to check for convergence, some solutions with a 
1 st order polynomial were pathological, and we adopt instead the 
stable 0 th order. 


For imaging, we first continuum subtracted the data by 
extracting the visibilities from the 210 MHz line free re¬ 
gion from 86.47 GHz to 86.67 GHz (rest frame) to define 
the continuum, channel averaging the visibilities, and 
then subtracting the channel averaged visibilities in the 
UV plane. The visibilities were imaged in MAPPING 
and the spectrum extracted in CASA (McMullin et al. 
2007) with an elliptical aperture 5.58" x 2.22", PA = 162° 
(three times the clean beam), centred on the phase cen¬ 
tre RA= 16 h 52 m 58 s .890, Dec= ±02°24'03".90. This in¬ 
cludes conversions from mJybeam -1 to mJy and allows 
for a slightly extended source. 

The noise in the imaged cube was estimated with the 
GILDAS task NOISE to be 0.2 mJy beam -1 . To estimate 
the noise in the spectrum we placed eight apertures, iden¬ 
tical to that used for the spectrum extraction, arranged 
around the science aperture so as to create a regular 3x3 
grid with the science spectrum in the centre. We then 
extracted the spectra from within these eight apertures. 
Since the cube is continuum subtracted, the signal in 
these apertures should be entirely noise (both thermal, 
and artificial from processing). For each spectrum the 
rms about zero was found, and then these eight rms’s 
were averaged to produce a final rms for the imaged spec¬ 
trum of 0.6 mJy channel -1 . The spectrum is shown in 
Figure 2. We fitted a single Gaussian to the SiO(2 — 1) 
line and find a reduced y 2 of 2.12. Attempts to fit a 
second Gaussian component consistently lead to patho¬ 
logical results due to all of the putative H 13 CN(1 — 0) 
flux lying in a single channel. 

The UV fitted and imaged spectra are compared in Fig¬ 
ure 3. The channel errors have been offset horizontally 
for clarity. Every channel is consistent between the two 
spectra within la. 

3.1. Continuum 

We used the 210 MHz region from 86.47 GHz to 
86.67 GHz (rest frame) to define the continuum. This 
region, expected to be line free, was verified as such by 
an inspection of the UV fitted spectrum. The imaged 
continuum is shown in Figure 4 and we find an aperture 
extracted flux density of 11.3 ± 1.7 mJy (using the same 
aperture and method as for the spectrum extraction). 
This is noticeably larger than the flux density of the el¬ 
liptical Gaussian fitted to the continuum UV visibilities, 
which gives a flux density of 9.70±0.29 mJy (9.7±1.5 mJy 
including the 15% absolute flux calibration uncertainty), 
although it is still consistent within ler. This is in part 
due to the extended irregularly shaped emission to the 
south-east of the pointing centre, which is not accounted 
for in the single component UV fit. Ideally we would fit 
further components to the residual visibilities (see e.g., 
Feruglio et al. 2013), but the poor signal to noise ra¬ 
tio of our data prevents meaningfully fitting additional 
components. 

An estimate of the continuum was also provided by 
the baseline fitted to the UV fitted spectrum. This mea¬ 
sured a continuum flux density of 9.58 ± 0.19 mJy, or, 
including our 15% absolute flux calibration uncertainty, 
9.6 ± 1.5 mJy. 

Our continuum flux density is slightly than ex¬ 
pected: Nakanishi et al. (2005) found 16.6±1.7mJy at 
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86.2 86.4 86.6 86.8 87.0 



Fig. 1. The UV fitted spectrum and the least squares two Gaussian plus baseline fit. Velocities are relative to the SiO(2 — 1) line centre, 
and rest frame frequencies in GHz are given at the top of the plot. The la channel uncertainties (median ±0.8 mJy) from the UV fitting 
are shown, and do not include the absolute flux calibration uncertainty. The horizontal bar indicates the region identified as line free and 
used for the continuum subtraction in the imaging step. The reduced \ 2 °f the fit is 1.52. 


Rest-frame frequency [GHz] 

86.2 86.4 86.6 86.8 87.0 



Fig. 2.— The imaged spectrum, extracted with an aperture as described in Section 3, and the least squares Gaussian fit to the SiO(2— 1) 
line. Velocities are relative to the SiO(2 — 1) line centre, and rest frame frequencies in GHz are given at the top of the plot. The la rms 
channel uncertainties (±0.6 mJy) are shown, and do not include the absolute flux calibration uncertainty. The horizontal bar indicates the 
region used for the continuum. The reduced x 2 °f the fit is 2.12. 
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Fig. 3.— A comparison of the UV fitted and Imaged spectra. Axes are as in Figures 1 and 2. The UV fitted spectrum has been continuum 
subtracted using the baseline fitted simultaneously with the Gaussian lines. Channel uncertainties have been offset left and right for clarity. 
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Fig. 4.— The NGC6240 3 mm continuum. Black crosses mark 
the positions of the AGNs as reported in Hagiwara et al. (2011). 
Contours are at -3, 3, 5, 10, 20 and 30(7, where a = 1.63 X 
10 -4 Jybeam -1 . 

~ 88.6 GHz rest frame and 10.8±1.1 mJy at ~ 110.2 GHz 
rest frame, with the Nobeyama Millimeter Array and 
RAINBOW interferometers. Their measurement is con¬ 
sistent with our imaged continuum at the 2 a level. Ex¬ 
panding the aperture used for the imaged continuum 
does not significantly change the recovered flux density. 


3.2. Spectral Lines 

Given the low signal to noise ratio of our observations we 
adopt an elementary statistical approach to line identifi¬ 
cation. A priori, we expect to detect H 13 CN(1 — 0) and 
H 13 CO + (l —0) based on previous detections of HCN and 
HCO + in NGC 6240, while we expect to detect SiO(2 —1) 
based on comparisons with other LIRGs and major merg¬ 
ers. We give parameters for these lines in Table 1. 

Only SiO(2 — 1) is clearly detected in both the UV 
fitted and imaged spectra. The H 13 CN(1 — 0) line ap¬ 
pears to be present in the UV fitted spectrum, but not 
clearly so in the imaged spectrum, and there is no sign 
of H 13 CO + (l — 0) in either. While, if detected, it would 
be at least partially blended with the SiO(2 — 1) line, 
there is absolutely no sign of emission redwards of the 
H 13 CO + (l — 0) line centre, or bluewards in the residuals 
of the SiO(2 — 1) line fit. 

For our statistical test of the line detections we treat 
the channel uncertainties as Gaussian. We sum the chan¬ 
nels within ±400 km s~ 1 of the line centre and sum the 
channel rms’s in quadrature, to create a single ‘detection’ 
channel and quantify this signal in units of ^detect) the 
rms of the detection channel. We define our upper limit 
on the flux density, as the flux density which has a 

2% chance of being measured at the value we see or less 
if the line is in fact present, i.e., 


Q 

'-'max 


^detect ± V^2cx^erf 1 (2p-l) 


( 7 ) 


channel, p is our 0.02 probability and erf is the Gauss 
error function. For p = 0.02 this gives S ma , x = detect + 
2.054a (Masci 2011) 8 . 


TABLE 1 

Observed Line Parameters 


Line 

Z'O 

[GHz] 

E u /k 

[K] 

Sline (UV) 
[Jykms 1 ] 

Sline (imaged) 
[Jy kms -1 ] 

H 13 CN(1 - 0) . a 

86.340 

4.14 

0.46 ± 0.07 

< 0.8 

H i3 CO+(l - 0) a 

86.754 

4.16 

< 0.4 

< 0.4 

SiO(2 - 1). b 

86.847 

6.25 

1.10 ±0.15 

0.90 ±0.09 


The uncertainties quoted here do not include the 15% absolute 
calibration uncertainty. 

“Spectral Line Atlas of Interstellar Molecules (SLAIM, Remijan 
et al. 2007). 

b The Cologne Database for Molecular Spectroscopy (CDMS, 
Muller et al. 2005). 


We use two measures to test whether the line is actually 
detected. Firstly, we test against the null hypothesis that 
we see a signal at least as bright as in the data, given 
that there is no line present - i.e., we assume that there 
is no chance of absorption. This is a one-tailed test with 
probability 


a = 


1 

2 


1 — erf 


x — p, 


( 8 ) 


and a value consistent with no line is 0.5. 

Secondly, we test against the null hypothesis that we 
see a signal with an absolute value at least as great as is 
seen, given that there is no line present - i.e, we allow for 
absorption. This is a two-tailed test with probability 


a = 1 — erf 


f \x- S | 

V V2a* 


(9) 


and a value consistent with no line is 1. While the one- 
tailed test specifically searches for emission lines, the two- 
tailed test makes no such distinction. 


For the aperture extracted SiO(2 — 1) line and contin¬ 
uum region we found one-tailed values of 7.0 x 10 -7 and 
0.53 respectively, precisely as expected for a clear line 
detection and a flat continuum. The two-tailed values 
are 1.4 x 10” 6 and 0.95, again consistent with a clear 
line detection and a flat continuum. 

For the imaged H 13 CN(1 — 0) line the one and two- 
tailed tests gave 0.014 and 0.028 respectively, while the 
UV fitted line gave 0.014 and 0.004. While not en¬ 
tirely conclusive, these tests argue for a detection of 
the H 13 CN(1 — 0) line in the UV fitted spectrum. For 
H 13 CO + (l —0) we subtracted the SiO(2 —1) line fit from 
the spectrum before running the tests. In all cases we 
find values completely consistent with no line detection 
with one and two-tailed tests giving 0.51 and 0.98 respec¬ 
tively. 

For the H 13 CN(1 — 0) line we find an upper limit 
S u dv < 0.8Jykms -1 from the imaged spectrum and 
^du < 0.9Jykins -1 from the UV fitting. The im¬ 
aged spectrum cannot be well fitted with a Gaussian 


where Sdetect is the measured signal in the detection 8 wise2.ipac.caltech.edu/staff/fmasci/UpperLimits_FM2011.pdf 
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due to almost all of the flux residing in one channel, 
but for the UV fitted spectrum the Gaussian fit to the 
line gives S v dv = 0.46 ± 0.07 Jykms -1 . The SiO(2 — 1) 
line was clearly detected in both spectra, and was fit¬ 
ted with line fluxes of S v dv = 1.10 ± 0.15 Jykms -1 and 
S v dv = 0.90 ± 0.09 Jykms -1 in the UV fitted and im¬ 
aged spectra respectively (not including the 15% abso¬ 
lute calibration uncertainty). For H 13 CO + (l — 0) we 
first subtracted the SiO(2 — 1) Gaussian line fit before 
summing the channel fluxes. We obtain an upper limit 
of S u dv < 0.4Jykms -1 for both spectra. Attempts 
to simultaneously fit two Gaussians to the SiO(2 — 1) 
and H 13 CO + (l — 0) line centres consistently reported 
zero line flux for H 13 CO + (l — 0), consistent with a non¬ 
detection/upper limit. The measured fluxes and upper 
limits are recorded in Table 1. 

We present 30MHz channels maps for H 13 CN(1 — 0), 
H 13 CO + (l —0) and SiO(2 —1) in Figure 5. H 13 CN(1 —0) 
is seen only in the channel containing the line centre, at 
a peak level of 3a. It is concentrated between the two 
nuclei, as was seen for HCN(1 — 0) by Nakanishi et al. 
(2005). There is no sign of H 13 CO + (l — 0) in any of the 
channels, consistent with the non-detection in the UV 
fitted and imaged spectra. Finally, SiO(2 — 1) appears 
at the 3 and 4cr level in the 0 and —110km s -1 channels 
centred on the southern nucleus, with a possible 3er de¬ 
tection in the +110 km s -1 channel about the northern 
nucleus. 

We briefly note than there is no sign of HN 13 C(1 — 0) in 
either of the spectra. Also, if the line were detected it 
would extend beyond the edge of our spectrum, hindering 
its analysis. We mark the line frequency in Figures 1 and 
2 for reference. 


(ALMA) HCN(4—3) observations of Scoville et al. (2015) 
have APEX single dish counterparts from Papadopoulos 
et al. (2014) and a flux ratio 0.90 ± 0.16, i.e., consistent 
with no spatial filtering/almost complete flux recovery. 
While extrapolating from the J = 4—3 transition to our 
J = 1 — 0 transitions is unacceptably unreliable, com¬ 
bined with the very compact HCN(1 — 0) emission and 
the comparisons in the previous paragraph this compar¬ 
ison does argue for our observations recovering most of 
the H 13 CN(1 - 0) flux. 

We note that when comparing ALMA observations 
with JCMT fluxes Scoville et al. (2015) raised signifi¬ 
cant concerns regarding the calibration accuracy of the 
single dish measurements. 

Wang et al. (2013) observed SiO(2 — 1) in a six lo¬ 
cal galaxies with the IRAM 30m, including NGC 6240. 
They found a double peaked line profile, with a com¬ 
bined flux of 0.72 ± 0.30 Jykms -19 , consistent with our 
observations. Wang et al. (2013) also observe what 
may be a marginal H 13 CO + detection, with a peak 
T* a ~ 0.55 mK. We approximate this as corresponding 
to cs 0.08 Jykms -1 , consistent with our upper limit on 
the line of 0.4 Jykms -1 . 

3.4. The LTE limit 

Papadopoulos et al. (2014) reported the discovery of an 
extremely high [ 12 C]/[ 13 C] abundance ratio in NGC 6240, 
evidenced by the [CO]/[ 13 CO] ratio of 300 - 500. Here we 
test whether this ratio is compatible with our observed 
H 13 CN(1 — 0) line flux under the assumption of local 
thermodynamic equilibrium (LTE). 

We compare our H 13 CN(1 —0) line flux to an aggregate 
of literature measurements of the HCN(1 — 0) line using 


3.3. Spatial Filtering and Missing Flux 

Interferometers sample discrete spatial frequencies, and 
as such may underestimate the true flux of a source. As 
there are no single dish observations of any of our tar¬ 
geted lines we have no direct means of measuring the frac¬ 
tion of recovered flux and instead use comparisons with 
other line and continuum measurements from the liter¬ 
ature. Nominally, our observations should only resolve 
out emission on scales larger than 12", three times larger 
than the largest line and continuum structures seen by 
Nakanishi et al. (2005), Iono et al. (2007) and Scoville 
et al. (2015), so it is unlikely that we have any signifi¬ 
cant loss of flux. However, given the detection of CO on 
scales ~ 25" by Feruglio et al. (2013) we proceed with 
caution, despite not expecting HCN or HCO + emission 
to be similarly extended. 

Nakanishi et al. (2005) found the HCN(1 —0) line to be 
more compact than both the continuum and HCO + (l — 
0). The HCO + (l — 0) line was slightly more extended, 
but still only out to scales ~ 4". It is of course possible 
that the observations of Nakanishi et al. (2005) are also 
filtering a significant fraction of the flux, and indeed their 
single dish flux to interferometric flux ratio was 0.68 for 
HCN(1 — 0), using the single dish data of Solomon et al. 
(1992). However, the single dish fluxes are extremely 
variable, and comparisons with more recent single dish 
observations give flux ratios of 1.1 ± 0.4 (Greve et al. 
2009) and 0.9 ± 0.4 (Krips et al. 2008). 

The high spatial resolution (0.5") interferometric 




T A dv 


2fcT A zA 
c 2 
he 3 


U=4 1 


8tt kv 2 


N„A 


ul 



1 - e -r 


r 


( 10 ) 

( 11 ) 


where £Ia,s is the solid angle of the primary beam and 
source respectively, N u is the column density of the 
species in level u in cm -2 and A u i is the Einstein A co¬ 
efficient for the transition between levels u and l (Gold¬ 
smith & Langer 1999). Combining these equations and 
assuming that HCN and H 13 CN share the same exci¬ 
tation temperature, T ex , so that Nhcn.j/Aj+scn j = 
[HCN]/[H 13 CN], we obtain for the J = 1 — 0 transition: 


Sidv! [HCN] -A 1:1 _ 0 T 2 (1 — e -Tl ) 

S 2 dv 2 [H 13 CN] A 2)1 _ 0 n (1 — e -T2 ) ’ 1 J 

where the subscripts 1 and 2 refer to HCN and H 13 CN re¬ 
spectively, Therefore, the lower limit on [HCN]/[H 13 CN] 
for a given flux ratio is found when both lines are opti¬ 
cally thin 10 . In this optically thin LTE case and using 


9 Converted from the TJ fluxes in their Table 1, using their 
quoted + e // = 95% and the 30m conversion factor S/T A = 
3.06+= ff/A e ff Jy K — , with A e ff = 0.63 from http://www. iram. 
es/IRAMES/mainWiki/Iram30mEfficiencies. 

10 Assuming that the optical depth of HCN(1 — 0) is always 
greater than that of H 13 CN(1 — 0) - a safe assumption. 










A a [arcsec] 


Fig. 5.— Continuum subtracted channel maps for H 13 CN(1 — 0), H 13 CO+(l — 0) and SiO(2 — 1). Black crosses mark the AGN as 
reported in Hagiwara et al. (2011) Contours are at —3, —2, 2, 3 and 4cr, where a = 3.9 x 10 —4 Jy beam -1 . H 13 CN(1 — 0) appears to be 
located between the two nuclei, co-located with the the HCN(l — 0) (Nakanishi et al. 2005). SiO(2 — 1) on the other hand appears to be 
concentrated around the more active southern nucleus. 
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the upper limit on our H 13 CN(1 — 0) line flux, we find 
[HCN]/[H 13 CN]> 12. However, exploring a range of op¬ 
tical depths we find that constrained by the J = 1 — 0 
lines alone the [HCN]/[H 13 CN] could be as large as 1000, 
consistent with the high values found by Papadopoulos 
et al. (2014). We further constrain the isotopologue ratio 
in Section 4 using LVG modelling. 

A high [ 12 C]/[ 13 C] ratio, leading to high [HCN]/[H 13 CN] 
and [HCO + ]/[H 13 CO + ] ratios, is also consistent with 
SiO(2 — 1) being brighter than both H 13 CN(1 — 0) 
and H 13 CO+(l - 0) (SiO(2 - l)/H 13 CO+(l - 0) > 3); 
SiO(2 — 1) is not over-luminous, rather the low abun¬ 
dances of the 13 C isotopologues leads to unusually faint 
H 13 CN(1 - 0) and H 13 CO+(l - 0) lines. Usero et al. 
(2004) found SiO(2 — l)/H 13 CO + (l — 0) intensity ratios 
of about 3 in some regions of NGC 1068, but this galaxy 
is unusually bright in SiO. 

4. LVG MODELLING 

We use the publicly available code RADEX 11 (van der 
Tak et al. 2007) with collisional data from the Leiden 
LAMBDA database (Green & Thaddeus 1974; Flower 
1999; Dayou & Balanga 2006; Lique et al. 2006; Du¬ 
mouchel et al. 2010; Yang et al. 2010) to model spectral 
line ratios as a function of the gas phase conditions. We 
employ an MCMC script using the Metropolis-Hastings 
algorithm (Metropolis et al. 1953; Hastings 1970) to ex¬ 
plore the parameter space. Details of our MCMC LVG 
model are included in Appendix A. We use this code to 
fit three gas phases and six species simultaneously. 

We combine our PdBI data with literature data for 
HCN and HCO + shown in Table 2 and adopt CO and 
13 CO data from Papadopoulos et al. (2014). We have 
available to us the complete HCN and HCO + ladders 
from J = 1 — 0toJ = 4 — 3 except for HCO + (2 — 1), 
the complete CO SLED up to and including J = 13 — 12, 
and the 13 CO J = 1 — 0, 2—1 and 6 — 5 lines. 

In collecting data from a range of telescopes and over a 
range of frequencies it is essential to consider whether 
there is a common putative source size and whether the 
single dish beam coupling to the source varies signifi¬ 
cantly between observations. We first convert the line 
fluxes from our observations and the literature into line 
luminosities: 


L'= 3.25xl0 7 2 fvr0 (1 + z )~ 3 ( J Av Svd * ) 

V GHz / \ Mpc J v Jy km s ) 


V Jy km s 1 
(13) 

where D l is the luminosity distance of the source and L' 
has units of Kkms” 1 pc 2 (Solomon et al. 1992). Ideally 
we would have a complete, spatially resolved sample of 
lines so that we could precisely compare the line bright¬ 
nesses at specific positions in the source, ensuring that 
we were comparing cospatial regions of equal size. Work¬ 
ing instead with almost entirely unresolved single dish 
observations we make the following approximations. 


Beam coupling: Interferometric observations of the 
dense gas tracers exist for only the J = 1 — 0 and 


TABLE 2 
LVG Line Inputs 


Line 

uo 

[GHz] 

•^line 

[Jy kms -1 ] 

References 3, 

HCN(1 - 0) 

88.63 

14 ±2 

1, 2, 3, 4 

HCN(2 - 1) 

177.3 

44 ± 7 

3 

HCN(3 — 2) 

265.9 

74 ± 7 

3, 4, 5 

HCN(4 — 3) 

354.5 

41 ±6 

5, 6 

H 13 CN(1 - 0) 

86.34 

0.46 ±0.07 

X 

HCO+(1-0) 

89.19 

21 ±3 

2 

H 13 CO+(l - 0) 

86.75 

< 0.44 

X 

HCO+(3-2) 

267.6 

141 ± 21 

5 

HCO+(4 — 3) 

356.7 

74 ± 9 

2, 5 


a l = Nakanishi et al. (2005), 2 = Greve et al. (2009), 3 = Krips 
et al. (2008), 4 = Gracia-Carpio et al. (2008), 5 = Papadopoulos 
et al. (2014), 6 = Scoville et al. (2015), x = this work. 

J = 4 — 3 HCN and HCO + lines. Convolved source 
sizes range from 4" to 1".5 as frequency increases, with 
the deconvolved low—J source sizes being quite uncer¬ 
tain, but with very compact high—J deconvolved sizes: 
l".l x 0".6 for HCN(4 — 3) (Scoville et al. 2015). Even if 
the low-J emission is more widespread it will still be much 
smaller than the single dish beam sizes at the associated 
sky frequencies, so that the beam coupling corrections 
are <1.05 and as such are much less significant than the 
single dish calibration uncertainties (> 10%). Also, this 
would be a systematic underestimate of the flux, so is 
accounted for to first order when using line ratios. 
Source size: In all of the observations, except for the 
HCN(4 — 3) of Scoville et al. (2015), the source size is 
so uncertain that we cannot make any convincing esti¬ 
mates of the source size beyond the flat assumption that 
there is a common size to all J levels. This could poten¬ 
tially introduce a bias into the models, with the lower—J 
lines being integrated over a larger area than the high—J. 
However, this acts in the opposite sense to any potential 
beam coupling errors, which reduce the measured low—J 
fluxes if the low—J lines are significantly more extended. 
To some extent these opposing errors correct for any in¬ 
accuracies introduced by assuming a common source size 
and not attempting to account for changes in the beam 
coupling. 

In addition to the observed line ratios we use the dynam¬ 
ical parameter K V i r to constrain our models (e.g., Greve 
et al. 2009; Papadopoulos et al. 2012, 2014), where 

_ du/dr _ 1.54 du / n H2 \”°- 5 , , 

vir (du/dr) vir ^ dr 11000 cm- 3 J ' 1 1 

The geometric coefficient a ranges from 1 — 2.5 with an 
expectation value {y/a) = 1.312. The velocity gradi¬ 
ent, du/dr, describes the change in line-of-sight velocity 
through the cloud in kms” 1 pc” 1 . In galactic GMCs in 
virial equilibrium K v j r ~ 1. However, as was pointed 
out by Papadopoulos et al. (2014), in the extreme en¬ 
vironments of LIRGs and ULIRGS a wider K vir range 
is necessary to capture the true environment due to un¬ 
usual GMCs and unbound gas. This is especially true 
given that we also model the shocked CO gas. We there¬ 
fore follow Papadopoulos et al. (2014) and adopt limits 
of 0.5 and 20.0 on the allowed K vlr values for our models. 


11 http://home.strw.leidenuniv.nl/~moldata/radex.html 


The modelling consisted of runs of increasing complex- 
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TABLE 3 

Three Phase Input Ranges. All values log 10 . 


Parameter 

n n 2 

dv/dr 

[CO]/[ 13 CO] 

[HC(N/Q +)1 

[h 13 c(n/o+j] 

^co 

^HCN 

^HCO+ 

/ 


Hot, shocked Diffuse Cold, dense 
2.0 - 3.5 1.0 - 3.0 0.5 - 1.5 

1-6 1-6 4-8 

-1-3- 

-1 - 3- 

1-3 

- 8 - -2 - 

-12-4 

-12-4 

- 2-0 - - 2-0 


ity to check for consistency both with existing data and 
between models. Here, we present only the final model. 
The precursors are included in the Appendix. 

The limits on the [Ai 2 c]/[X 13 c] range were taken from 
the analysis of the HCN(1 —0)/H 13 CN(l —0) line ratio in 
Section 3.4 to be 10 and 1000, while the limits on Tk, «h 2 
and du/d r are chosen to be consistent with Papadopoulos 
et al. (2014), although we extend the lower limit on du/d r 
to 0.1kms _1 pc _1 to better accommodate diffuse CO. 
The temperature of the background blackbody, T^ g , is 
set to 3K. For models with free abundances we adopted 
—12 < log 10 (X mo i) < —4 for HCN and HCO + and —8 < 
log 10 (*co) < -2 for CO. 


there is a free [CO]/[ 13 CO] ratio. In the cold, dense 
phase there is also an [HCN]/[H 13 CN] ratio, which is 
shared with HCO + . Inspired by Zhang et al. (2014) and 
Kamenetzky et al. (2014) we relate the fluxes from the 
three phases through two free scaling factors /1 and / 2 , 
such that: 


^ ^shocked fl T ^diffuse + ^ dense f 2 • (15) 

It is important to appreciate that the /’s are NOT 
the fractional flux densities emerging from the phases, 
as the fluxes from each phase may intrinsically be much 
greater (or smaller) than Sdiffuse- They are merely free 
parameters to allow the MCMC to scale the significance 
of the two phases. Physically, they are indirectly related 
to the beam filling factors of the phases, whereas the 
formulation of Kamenetzky et al. (2014), who fitted line 
intensities instead of line ratios, used two beam filling 
factors directly as free parameters. 

The contrast factors allow us to abstract away the issue 
that different molecular species present different beam 
filling factors. For the unresolved single dish observa¬ 
tions this is essential, but even for high resolution in¬ 
terferometric observations this is a desirable property, 
allowing for that fact that the ISM can vary significantly 
over parsec scales (barely resolvable with ALMA in even 
the closest galaxies). 


4.1. Three Phase Modelling 

A concern with simultaneously modelling CO, HCN and 
HCO + is that about 30% of the CO(l — 0) line flux is far 
more extended that any other molecular species, tracing 
an outflow out to almost 30", while the HCN and HCO + 
emission only extends out to < 2". (Feruglio et al. 2013; 
Nakanishi et al. 2005), potentially leading to significant 
errors and bias if we force it into the same gas phase 
as HCN and HCO + . We therefore adopt a three phase 
model consisting of a diffuse, extended CO phase, a hot 
and shocked CO phase and a cold, dense multi-species 
phase. 

The parameter ranges are given in Table 3. The tem¬ 
perature and density ranges are chosen so as to force the 
phases to be hot, diffuse or dense respectively, while still 
allowing a large range for each parameter to explore. 

HCN and HCO + are only included in the cold, dense 
phase. This is motivated by the absence of any HCN 
emission from the shocked phase in the two phase mod¬ 
els (Appendix C) . There is a potential issue with exclud¬ 
ing HCO + from the shocked and/or diffuse phase as it 
shows evidence in the precursor models of having a small 
but significant contribution from hotter and more diffuse 
gas. The choice to restrict HCO + to the dense phase is a 
pragmatic one: we simply do not have enough observed 
lines to constrain a more accurate model of HCO + . This 
should not significantly affect the results as the CO lines 
are the primary determinants of the conditions in the dif¬ 
fuse and shocked phases, and the warmer phase HCO + 
emission is less than 10% of the flux in any of the lines. 
Furthermore, it is very unlikely to be present in signifi¬ 
cant quantities in the hot, shocked gas, nor in the fully 
extended diffuse phase; rather it is expected on the UV 
illuminated edges of molecular clouds, which we do not 
model at all. 

CO is present in all three phases, and in each phase 


The results of the three phase modelling are shown in 
Figures 6, 7 and 8, while the numerical values are tabu¬ 
lated in Table 4. The model successfully splits into three 
clear phases, without any tension against the parame¬ 
ter ranges. Xco varies greatly between the gas phases, 
being greatest in the diffuse phase lying near the canoni¬ 
cal value of 2 x 10 -4 , before dropping slightly to around 
1 x 10 -5 in the shocked phase and then dropping rapidly 
to around 10~ 6 — 10~ 7 in the cold, dense phase; strongly 
suggestive of freeze-out of CO onto dust grains in the 
cold, dense phase. 

The isotopologue abundance ratios are also largely as 
we would expect: both the shocked and the dense gas 
phases have relatively low [CO]/[ 13 CO] ratios of 30—100, 
likely representative of the true [ 12 C]/[ 13 C] ratio in the 
galaxy. The diffuse phase present a slightly elevated ratio 
of ~ 100 — 200, consistent with this ubiquitous gas phase 
experiencing selective photodissociation due to the nu¬ 
clear starbursts. Finally, the [HCN]/[H 13 CN] ratio (re¬ 
call that this is shared with HCO + ) is the highest, at 
around 300, consistent with isotope fractionation in cold, 
dense molecular cores. 

The SLEDs in Figure 7 show exceptionally good fits to 
the observed lines. The cold, dense phase has very little 
contribution to the CO or 13 CO SLEDs, and in both of 
these the shocked phase dominates the emission at J > 4, 
while the low—J lines are set by the ubiquitous diffuse 
phase. 

For each point in the MCMC trace we can calculate 
K v i r : the K V i r pdfs are shown in Figure 8. Both the 
diffuse and dense phases present doubly peaked pdfs, in 
particular with a peak about K vir = 1, while the shocked 
phase pdf steadily increases towards the upper limit of 
K vir = 20. K v i r is generally a poorly constrained pa¬ 
rameter in our models, but these results are suggestive of 
approximately virialised gas dominating the diffuse and 
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TABLE 4 

Three Phase Results. All values log 10 . Values in 

PARENTHESIS ARE THE BEST FIT VALUES. 


Parameter 

Hot, shocked 

Diffuse 

Cold, dense 

T k 

3.2g; a (3.2) 

2.3 2 ;§(2.4) 

i-il:g(o.9) 

n H 2 

3.6|;1(3.5) 

2-51:5(2.4) 

5.61:3(6.6) 

dv/dr 

0.8,!,;° (0.2) 

9-1— 0.6 (9-2) 

1-8} : ?(1.5) 

[CO] 

pcoJ 

1.6j;§(1.8) 

2.2 2 ;|(2.2) 

2.2 2 ;1(2.5) 

[HC(N/0+)j 

[H 1 3C(N/0+)] 

- 

- 

2-51:1(2.6) 

- Y CO 

-5,o: 4 ;®(-5.i) 

-3.8I 2 ;I(-3.7) 

-6.6I®;«(-7.4) 

-1'hcn 

- 

- 

-8.0l£j(-9.6) 

-Wco+ 

- 

- 

— 8-8=i;g(—10.2) 

/ 

—0.5z8;|(— 0 . 6 ) 

- 

—0-7:?;g(—0.2) 

Mass [M 0 ] a 

9-4l:|(9.1) 

9-31:1(8.7) 

10-3g°6 8 (10.0) 


a CO estimated mass from calculated aco and the observed 
CO(l — 0) line. 


dense gas phases, while the shocked phase is highly tur¬ 
bulent. 

We also show the [HCN]/[HCO + ] abundance ratio pdf 
in Figure 9. This ratio is very tightly constrained as 
6 - 21 : 1 - 


We can also calculate a mo i(i-o)i where a mo i(i_o) = 
Mdense/^moi(i-o) the ra -ti° of the gas mass of the phase 
to the J = 1 — 0 line luminosity of the species from that 
phase 12 . We use equation A4 for optically thick gas from 
Papadopoulos et al. (2012) 13 to derive a mo i(i-o) f° r CO 
and HCN and present the pdfs, along with the corre¬ 
sponding masses, in Figure 10. The spread in each of 
the aco pdfs is between 1.5 and 2 dex, with the centres 
shifting to increasing aco as we move from the diffuse 
phase, through the shocked phase and onto the dense 
phase. There is less spread in the masses, and in partic¬ 
ular the HCN derived and CO derived dense gas masses 
are very consistent with almost identical pdfs but for a 
systematic shift of ~ —0.2 dex on the CO mass estimates. 
These are discussed in Section 5.3. 

Interestingly, the thermal pressure of the shocked and 
dense gas phases is P/k ~ 10' Kkms -1 - the same value 
that Dopita et al. (2005) found necessary to fit the SED 
in NGC6240 and the value found for the hot, diffuse 
phase in Kamenetzky et al. (2014). This makes sense 
in light of the three phase model: the diffuse CO and 
shocked CO dominate the CO SLED, so when Kamenet¬ 
zky et al. (2014) fitted to the CO SLED only they recov¬ 
ered these phases. However, the diffuse phase is only 
partially co-located with the shocked phase so is not 
in pressure balance, whereas the dense phase, identified 
here by fitting HCN, HCO + and CO simultaneously, is 
co-located with the shocked phase, and approximately in 
pressure balance with it. On the other hand, we would 
expect a significant fraction of the pressure balance to 
be due to macroscopic turbulence, not just thermal pres¬ 
sure. 
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4.2. Comparison to the Literature 

The dense, HCN bearing molecular gas in NGC 6240 has 
been modelled previously, by Krips et al. (2008), Greve 
et al. (2009) and Papadopoulos et al. (2014) (Kamenet¬ 
zky et al. 2014 only modelled the CO SLED). While 
Krips et al. (2008) modelled a single gas phase for both 
HCN and HCO + , Greve et al. (2009) and Papadopoulos 
et al. (2014) fitted HCN and HCO + independently, each 
with their own gas phase parameters. 

While Krips et al. (2008) and Greve et al. (2009) both 
found similar gas kinetic temperature ranges (20 — 120 K 
and 60 — 120 K respctively), Greve et al. (2009) found 
much larger tth 2 in the HCN phase of ~ 10 5 cm -3 , cf. 
10 3 - 5 - 10 4 ' 5 cm -3 found by Krips et al. (2008). On the 
other hand, the HCO + phase had an ?lh 2 of ~ 10 4 cm -3 , 
consistent with the results of Krips et al. (2008). 

Papadopoulos et al. (2014) present their results as 2D 
pdfs, finding degenerate “bananas” ranging from «h 2 — 
10 6 cm- 3 /T k ~ 8 K to n fl2 ~ 10 4 cm- 3 /T k ~ 1000 K, al¬ 
though the HCN pdf appears to peak at the high-density, 
low-T k end of the “banana”. 

Our results are broadly consistent with these previ¬ 
ous works. Both our precursor models and the full three 
phase model argue for relatively cool (8 — 30 K), dense 
( 10 4 - 5 _ 10 6 - 5 cm -3 ) gas dominating the HCN and HCO + 
emission. These models echo the trend seen by Greve 
et al. (2009) and Papadopoulos et al. (2014) of HCO + 
preferring a ~ 0.5 dex lower nn 2 than HCN while shar¬ 
ing a similar kinetic temperature. Our [HCN]/[HCO + ] 
abundance ratio of 5 . 3,4 7 is consistent with Krips et al. 
(2008) who found “around 10 ”. 

The multiple phases of the molecular gas in NGC 6240 
have been characterised previously, including by Ar- 
mus et al. (2006) with Spitzer observations of H 2 vibra¬ 
tional lines, Papadopoulos et al. (2014) with the afore¬ 
mentioned spectral decomposition of the CO SLED and 
Kamenetzky et al. (2014) with a two phase MCMC model 
of the CO SLED. Our models are the first to simulta¬ 
neously fit multiple species in multiple gas phases. As 
was pointed out by Kamenetzky et al. (2014), and as 
we found here comparing to the results of Papadopou¬ 
los et al. (2014), it is important to simultaneously fit 
multiple components to SLEDs to prevent the biasing of 
results. 

From modelling of H 2 vibrational lines Armus et al. 
(2006) found two highly excited phases of the molecular 
gas, with 6.7 x 10 6 Mq at T k = 957 K and 1.6 x 10 9 Mq 
at T k = 164 K respectively. Papadopoulos et al. (2014) 
first fitted the HCN and HCO + SLEDs, before using 
these results to fit high—J components of the CO SLED, 
then fitting the residual CO SLED. This produces three 
gas phases with densities log 10 (uh 2 ) = 5.0, 4.3 and 3.0 
and kinetic temperatures log 10 (T k ) = 1.5, 2.6 and 2.0. 
The model of Kamenetzky et al. (2014) used a nested 
sampling algorithm to fit two gas phases to the CO 
SLED, finding a hot and a cold phase with best fit den¬ 
sities log 10 (nn 2 ) = 4.1 and 5.8 and kinetic temperatures 
log 10 (T k ) = 3.1 and 1.2. These are reassuringly similar 
to the hot, shocked and cold, dense phase results of our 
multi-species, multi-phase analysis. It is curious however 
that the cold phase of their two phase model was closer 
to the fainter cold, dense phase than the brighter diffuse 
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Fig. 6. — The step plot for the three phase HCN+HCO++CO model results. Contours are at the 68%, 95% and 99% credible intervals. 
Red, black and blue contours are the shocked, diffuse and dense gas phases respectively. In the \Xi 2 q\/[Xi 3 c ] plots the green contours 
correspond to the [HCN]/[H 13 CN] abundance ratio (shared with HCO"*"), and the blue the [CO]/[ 13 CO] ratio, both in the cold, dense 
phase. 


CO phase. There is no clear explanation for this, but 
it could very possibly be an artefact of the change from 
two to three phases. 

5. DISCUSSION 

5.1. Validity of Multi-species Models 

An important question of our analysis is whether or not 
it is appropriate to model multiple species as sharing 
the same gas phase. In reality there is a continuum 
of conditions within any telescope beam and any phase 
based model is only an approximation of the conditions 


dominating the emission, which will vary from molecule 
to molecule. Nevertheless, combining molecular species 
which are expected to trace similar regions appears to 
be a powerful tool for estimating molecular abundance 
ratios and gas conditions. 

The combination of two or more species in a single 
model is a compromise, with the extent of the assump¬ 
tions dependent upon how similar the regions traced 
by the different species are likely to be. In our three 
phase model, the dense phase is shared by HCN, HCO + 
and CO, as well as their 13 C isotopologues. That CO 
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Fig. 7.— SLED fits for the three phase HCN+HCO + +CO model, showing the relative contribution of each gas phase. The same 
normalisation has been used for the CO and 13 CO SLEDs. 


is present alongside other molecules goes without ques¬ 
tion, even if its abundance is reduced due to freeze-out 
onto ice grains, and the other two gas phases prevent CO 
from being heavily biased by the dense conditions; the 
colocation of HCN and HCO + is less concrete. While 
HCN and HCO + are almost certainly colocated, HCO + 
is quite possibly more extended, tracing warmer and less 
UV shielded regions (Fuente et al. 2005). On the other 
hand, it is not expected to be found in the very diffuse 
and extended CO gas phase we use here. Unfortunately, 
while an additional phase would be desirable we would 
not be able to place meaningful constraints upon the ad¬ 
ditional phase due to the paucity of lines. Two phase 
HCN and HCO + models (Appendix C.3) suggest that the 
contribution of HCO + in a slightly hotter and less dense 


phase (physically this would most likely correspond to 
the transition region between the hot, diffuse phase and 
the cold, dense cores) is less than 10% of any HCO + 
line. We believe that the colocation of HCN and HCO + 
in the cold, dense phase is therefore an undesirable but 
necessary, and acceptable, compromise. 

5.2. The / 12 Cy// 13 Cy Abundance Ratio in NGC6240 

Papadopoulos et al. (2014) constrained the high—J CO 
lines using the LVG solutions for HCN and HCO + before 
fitting a third phase to the remaining CO lines. When 
they then applied these three phases to their 13 CO data, 
a very high [CO]/[ 13 CO] abundance ratio of 300 — 500 
was required for consistency and they argued for that 
this was evidence for a similarly elevated [ 12 C]/[ 13 C] ra- 
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Fig. 8. — The K v i r pdfs for the three phase model. Colours are 
as in Figure 6, with red, black and blue correspond to the shocked, 
diffuse and dense phases respectively. Both the diffuse and dense 
phases present doubly peaked pdfs, in particular with a peak about 
K v i r = 1, while the shocked phase pdf steadily increases towards 
the upper limit of K v i r = 20. 



!og 10 ([HCN]/[HCO + ]) 

Fig. 9.— The [HCN]/[HCO + ] abundance ratio pdf from the three 
phase model. The value is tightly constrained as 6.2p|. 

tio. Here, we have used our new observations of H 13 CN 
and H 13 CO + in combination with the CO and 13 CO 
lines from the literature to place new constraints on the 
[ 12 C]/[ 13 C] abundance ratio. We have found that instead 
of arising in the HCN and HCO + bearing dense gas, the 
high—J CO lines are due to the hot, shocked gas phase, 
which in turn changes the implied [CO]/[ 13 CO] ratios. 

A naive approach, simply taking the 68% credible in¬ 
tervals on [CO]/[ 13 CO] and [HCN]/[H 13 CN] as the brack¬ 
ets on the [ 12 C]/[ 13 C] ratio, places the ratio somewhere 
between 100 and 1000: however, this does not account 
for the much greater abundance of CO wrt HCN nor the 
dominance by mass of the dense phase. 

We adopt a mass and abundance weighted average and 
find [ 12 C]/[ 13 C] = 986 30 . The pdf of the ratio is shown 
in Figure 11. The ratio is strongly peaked at “normal” 
ULIRG values of about 100, although there is an ex¬ 
tended tail out to 1000. The model is compatible within 
— lcr with a Galactic value, while the +3er boundary is at 
550. We therefore argue against a super-ULIRG elevated 
[ 12 C]/[ 13 C] abundance ratio, but cannot conclusively ex¬ 
clude the range found by Papadopoulos et al. (2014). 
The ideal solution, a full chemical-liydrodynamical model 


with radiation transfer, is beyond the scope of this pa¬ 
per. Nevertheless, it appears that the [ 12 C]/[ 13 C] ratio 
in NGC6240 is most similar to other ULIRGs. 

Are the derived isotopologue abundance ratios consis¬ 
tent with the chemical models presented in Section 1.1? 
Given the wide range of results from the chemical models, 
and their strong dependences on temperature, density, 
metallicity, the C/O ratio and the H 2 ortho-para ratio it 
is hard to say with certainty. We find [HCN]/[H 13 CN] 
~ 2x [CO]/[ 13 CO] in the dense and diffuse gas phases, 
consistent with Roueff et al. (2015). The rather high 
[HCO + ]/[H 13 CO + ] (which in the three phase model is 
fixed to be the same as [HCN]/[H 13 CN]) is surprising, 
but not so extreme as to cause concern; Langer et al. 
(1984) found [HCO+]/[H 13 CO+] ~ 1.5x [CO]/[ 13 CO] in 
10K, 10 4 cm -3 , low metallicity gas. Furthermore, since 
this ratio is fixed to [HCN]/[H 13 CN] this could be artifi¬ 
cially elevating it (see Figure 17 in the Appendix). 

The [CO]/[ 13 CO] ratio appears to increase from 40|{j in 
the hot phase to IGO^q 0 in the cold, dense phase. There 
is slight tension between these results and our interpreta¬ 
tion of the model results as being evidence for ICE: in the 
case of ICE fractionation in the cold, dense phase should 
lower the [CO]/[ 13 CO] ratio; it is therefore contradictory 
if the hot, shocked phase, where there is no isotope frac¬ 
tionation, presents a lower [CO]/[ 13 CO] ratio than the 
dense phase. However, there is still a region of parame¬ 
ter space within the 68% credible interval about the mean 
values where the [CO]/[ 13 CO] ratio is lower in the dense 
phase, so while there is tension it is not that great, and 
the model is still consistent with isotope fractionation. 
Furthermore, the hot phase [CO]/[ 13 CO] is almost en¬ 
tirely determined by the highly uncertain 13 CO J = 6 — 5 
line, while the dense phase CO and 13 CO emission is al¬ 
most zero (Figure 7), so that the [CO]/[ 13 CO] ratio is 
highly susceptible to small errors. The [CO]/[ 13 CO] vs T k 
panel of Figure 6 shows both that there is a large spread 
in all of the [CO]/[ 13 CO] values, and that there is a very 
large overlapping region within the 68% credible inter¬ 
vals. Interestingly, within this spread in [CO]/[ 13 CO] all 
three phases are consistent with a value of ~ 80, closer 
to the ~ 70 of the Milky Way than to a globally elevated 
elemental [ 12 C]/[ 13 C] ratio. 

It is possible that due to insufficient UV-plane sampling 
we have lost some of the continuum and line fluxes, as 
was mentioned in Section 3.1. The effect of this would 
be to bias our [HCN]/[H 13 CN] and [HCO+]/[H 13 CO+] to 
higher values; our finding that the elemental [ 12 C]/[ 13 C] 
ratio is not globally elevated in NGC6240 is therefore 
robust against this potential source of error, as the error 
would push the results towards an elevated value. 

5.3. Line Luminosity to Gas Mass Conversion Factors 

The ratio between line luminosity and gas mass, a , is of 
great interest as it provides a means of estimating galaxy 
gas fractions from observations of a single molecular line, 
usually CO. However, the determination of a general cteo 
is particularly complicated, with dependences on lumi¬ 
nosity (Solomon et al. 1987) and metallicity, leading to 
a range of values from 0.3 to almost 300, and a typical 
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Fig. 10. — The ago factor for each of the gas phases of the three phase model and the «hcn (encompassing the dense phase) (left), and 
the corresponding mass of the gas phase (right). The diffuse and shocked phases are consistent with the “standard” oiqq values of between 
0.5 and 4, while the dense gas phase is considerably higher. While contributing a very small fraction of the CO(l — 0) line flux the dense 
gas phase dominates the molecular gas mass in the galaxy. The CO estimate and the HCN estimate of the mass of the dense phase are 
highly consistent. 



Fig. 11. - The pdf for the [ 12 C]/[ 13 C] ratio in NGC 6240, derived 
from the three phase model in Section C.5 using a mass and abun¬ 
dance weighted average. The ratio shows signs of being slightly 
elevated, but the upper 68% credible interval is less than half the 
value found by Papadopoulos et al. (2014) and the pdf is most 
consistent with the “normal” ULIRG abundance ratio. 

Milky Way value of 4.8 dropping to about 1 in z = 2 
SMGs and 0.8 in (U)LIRGs 14 (Downes & Solomon 1998; 
Bolatto et al. 2013). Nevertheless, it has recently been 
argued that aco may be underestimated in (U)LIRGs, 
and a Galactic value may be more accurate (Papadopou¬ 
los et al. 2012), due to gas over-densities arising from 
over-pressurisation due to turbulence. 

There is interest in «hcn as a possible estimator of 
the dense molecular gas mass, which may be more ro- 

14 Downes & Solomon (1998) found a range for aco from 0.3 to 
1.3 in ULIRGs. 


bust over a wide range of luminosities and may more 
accurately trace the specifically star forming gas (Gao 
& Solomon 2004; Greve et al. 2009; Garcia-Burillo et al. 
2012). This in turn has been complicated by high resolu¬ 
tion studies of ULIRGs which suggest that HCN emission 
may be dominated by shocked, X-ray heated or infra¬ 
red pumped emission (Aalto et al. 2007; Sakamoto et al. 
2010; Aalto et al. 2012; Garcia-Burillo et al. 2014; Viti 
et al. 2014; Aalto et al. 2015; Martin et al. 2015; Tunnard 
et al. 2015). 

We calculate a mo i and corresponding mass pdfs for 
each gas phase. The values for our three phase model are 
shown in Figure 10 where we plot aco for each of the 
three gas phases, and «hcn for the dense phase, as well 
as the corresponding masses. The aco pdfs are broad, 
with spans of about 1.5 — 2 dex. The diffuse and shocked 
phases both cover the range aco = 1 — 10, while the dif¬ 
fuse phase extends down to 0.1. The dense phase aco is 
much higher, due to a lower A'co> higher u.h 2 and higher 
du/d r in this phase, all of which contribute to an in¬ 
creased aco- Despite contributing < 5% of the CO (1 — 0) 
line flux, the dense phase accounts for almost all of the 
molecular gas mass. The dense phase gas masses derived 
using HCN and using CO are consistent but for an un¬ 
explained systematic offset of —0.2 dex on the CO mass 
estimate. Ohcn = 32f§: higher than the canonical value 
of 10 but similar to the 13 — 31.1 found by Papadopoulos 
et al. (2014), although some regions of their Figure 14 
do extend to ohcn = 60. Gao & Solomon (2004) pre¬ 
dicted ohcn ~ 25 for Ahcn = 2 x 10 -8 and du/d r = 
5kms _1 pc -1 : our lower Ahcn and higher du/d r are re¬ 
sponsible for our slightly higher ancN- While the + la- 
level is very high, the best fit and —la level are consis¬ 
tent with Greve et al. (2009), who found ancN = 17 — 37 
in NGC 6240. 

We derive a global aco by calculating the sum of the 
masses of each gas phase, and dividing by the total 



























16 



lo glo(«CO) 


Fig. 12.— The global aco for NGC6240 derived from the three 
phase model with the corresponding total molecular gas mass. 


CO (1 — 0) line luminosity. The resultant aco pdf and the 
corresponding mass pdf for NGC 6240 is shown in Figure 
12. This global aco is much more tightly constrained 
than the those of the individual phases, and we find 
a C o = 1-5H (l°gio (M/[Mq]) = lO.llg g): consistent 
with a Galactic aco, as predicted by Papadopoulos et al. 
( 2012 ), but perhaps more consistent with the canonical 
aco of starbursts (0.8, Downes & Solomon 1998). No¬ 
tably, the very low (CO) luminosity dense phase con¬ 
tributes the majority of the gas mass. 

The high ancN we derive is contrary to the findings of 
Garcia-Burillo et al. (2012), who argued for a lower Ohcn 
in (U)LIRGs, matching the lower aco- It is very possible 
that this contradiction is simply due to the extremely 
unusual properties of NGC 6240; NGC 6240 is perhaps 
not even representative of other (U)LIRGs, existing in a 
class of its own. 


column densities of the clouds are ~ 10 23 cm -2 then we 
are instead seeing multiple small ( 0.1 — 0.01 pc) dense 
clouds. 

We therefore propose the following situation. Within 
the turbulent internuclear disk instabilities lead to rapid, 
localised cooling of pockets of gas, leading to multiple 
small, cold, dense clouds embedded in pressure equilib¬ 
rium with the hot, diffuse phase but radiatively decou¬ 
pled due to turbulent Doppler shifting of the lines. This 
rapid cooling leads to an evolving chemistry, potentially 
with molecules less stable than CO destroyed in the high 
temperatures of the shocks reforming as the gas cools. 
These clouds are still turbulent however, and may indeed 
be transient, explaining the lack of star formation. 

In the three phase model the CO abundance is signif¬ 
icantly lower in the cold, dense phase than in the hot 
phase or diffuse phase: we have almost a continuum of 
Xco from the the diffuse phase value of ~ 10 ~ 4 , falling 
to ~ 10 -5 in the shocked phase and then to ~ ICC 6 
in the dense phase. This is consistent with the many 
cold clouds interpretation above, with freeze out of CO 
onto dust grain mantles depleting the gas phase CO in 
cold, dense clouds - as observed in Milky Way cold cores 
(Kramer et al. 1999; Tafalla et al. 2002; Liu et al. 2013; 
Ripple et al. 2013). This might also explain the very 
strangely hot and diffuse CS model solutions from the 
precursor models (Table 7): CS is depleted along with 
CO (Tafalla et al. 2002), so the only CS we observe is 
in the hot gas phase. Unfortunately this is highly spec¬ 
ulative, and cannot be confirmed without much higher 
spatial resolution observations able to resolve the high 
density molecular clouds. 

We note that with njj 2 ~ 1 x 10 5 , a uniform density 
profile and spherical gas distribution the dynamical mass 
implies an upper limit on the dense gas volume filling 
factor of 0.002 within the central 600 pc of NGC 6240. 
This is consistent with the image, presented above, of 
the molecular disk being almost entirely filled with hot, 
shocked gas, but with the majority (by mass) of the 
molecular gas surviving in cold, dense clouds embedded 
within the hot phase. 


5.4. Dense Gas and Star Formation in NGC6240 

The results of our LVG models present an intriguing pic¬ 
ture of the dense molecular gas in NGC 6240. We appear 
to be finding large quantities of very cold (~ 8 K) gas 
within the internuclear molecular disk, which is also host 
to extensively shocked and very hot (~ 2000 K) molecu¬ 
lar gas. The question is how can such a large quantity of 
cold, dense gas coexist with the hot phase. Furthermore, 
star formation in NGC 6240 is not found in the internu¬ 
clear molecular disk; it is instead concentrated in the two 
nuclei. If we truly have almost 10 10 M 0 of cold, dense 
gas then we must be able to explain why this gas is not 
actively forming stars. 

An insight is offered if we apply an upper limit to the 
H 2 column density of the dense phase. There is an inher¬ 
ent degeneracy within the LVG model between a single, 
large cloud and multiple, small clouds distributed in ve¬ 
locity space across the large scale velocity field of the 
galaxy. However, for the dense gas phase a single cloud 
implies column densities in excess of 10 26 cm -2 , which 
are extremely unlikely. If instead we assume that the 


Since the majority of the dense gas and the ongoing star 
formation are separate in NGC 6240 (Tacconi et al. 1999) 
we cannot calculate reliable local gas depletion times 
and instead we explore the galaxy averaged gas deple¬ 
tion time. Combining the star formation rate (SFR) of 
60±30 M 0 yr _1 (Yun & Carilli 2002; Feruglio et al. 2013) 
with our Ohcn(i-o) an d -^hcn(i-o) we can °kt a i n es ti~ 
mates of the gas depletion time in NGC 6240, and we 
find a depletion time 64 s x 10 s yr. 

However, the dynamical mass in the central 600 pc of 
~ 10 10 M 0 is significantly less than the dense gas mass 
implied by HCN and the dense phase CO. If we exclude 
results predicting gas masses in excess of the dynamical 
mass the depletion time is reduced to 0.9j;g x 10 8 yr. 

We do however question whether the dynamical mass 
is in fact particularly meaningful in the central regions of 
NGC 6240. This galaxy is an ongoing major merger, so 
it is not clear that the velocity dispersion should actually 
correlate with the mass. Estimating the mass from the 
velocity dispersion assumes either virialisation or stable 
orbital rotation, neither of which is necessarily true for 
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the central gas in NGC6240, which is dissipating turbu¬ 
lent energy from the merger (Tacconi et al. 1999) and is 
most likely a transient feature itself. It is very possible 
that this explains the super-dynamical masses suggested 
by some regions of our results and a large region of the 
parameter space in Papadopoulos et al. (2014). 

5.5. Implications of the [HCN]/[HCO + ] Abundance 

Ratio 

We have found [HCN]/[HCO+] = 6.2j;jj in NGC6240. 
This compares favourably to the value of “around 10” 
found by Krips et al. (2008). NGC6240 is frequently 
compared to Arp 220, where Tunnard et al. (2015) found 
[HCN]/[HCO + > 18 in the western nucleus (WN) and 
[HCN /[HCO + 0.5 — 5 in the eastern nucleus (EN). The 
abundance ratio for NGC 6240 appears to be most con¬ 
sistent with the less luminous EN, which is dominated 
by star formation (whereas Arp 220 WN may contain an 
energetically significant AGN). 

However, in NGC 6240 the comparison is probably less 
meaningful. Unlike in most of the galaxies of Krips et al. 
(2008), the majority of the molecular gas is spatially sep¬ 
arated from the starbursts and putative AGNs in the two 
nuclei. As Krips et al. (2008) pointed out, the extreme 
conditions in NGC 6240 make it quite unique, and the 
abundance ratios may not carry the same meaning as 
they do in other galaxies and Arp 220 (which may it¬ 
self be in a subset of galaxies possessing Compact Ob¬ 
scured Nuclei (CONs), which generate unique chemistries 
Costagliola et al. 2015). In the context of Krips et al. 
(2008), our results place NGC 6240 in the intermediate 
region between a starburst and AGN dominated galaxy. 
The derived n^ 2 is higher than the usual < 10 4 5 cm' 3 
for AGNs, while the [HCN]/[HCO + ] ratio is higher than 
the 0.01 — 1 of starburst but lower than the > 10 of 
AGN. Given the separation of the starburst and AGN 
nuclei from the bulk of the molecular gas (which lies be¬ 
tween the two nuclei) the fact that NGC 6240 does not 
fit clearly into any specific category is not surprising. 

6. CONCLUSIONS 

We have used new observations of the 13 C isotopologue 
lines of HCN and HCO + in NGC 6240 to model the physi¬ 
cal conditions and abundance ratios of the molecular gas. 
We have: 

• presented the first observations of SiO(2— 1) and of 
the isotopologue line H 13 CN(1 — 0) in NGC 6240, 
and an upper limit on H 13 CO + (l — 0). 

• combined these observations with literature data 
as inputs for an MCMC wrapper for RADEX, 
modelling three gas phases with six species simul¬ 
taneously to account for shocked molecular gas, 
extended diffuse molecular gas and cold, dense 
molecular gas, to find the conditions of the gas in 
NGC 6240. These models suggest that the cold, 
dense gas exists in multiple, small (0.01-0.1 pc) 
clouds embedded in the hot, shocked inter-nuclear 
gas disk. The cold gas is particularly cold with 
Tk = 13 2 q, while the hot phase reaches tempera¬ 
tures of almost 2000 K. 

• used our MCMC code and the new 13 C observa¬ 
tions to demonstrate that the very high [ 12 C]/[ 13 C] 


ratio in NGC 6240 found by Papadopoulos et al. 
(2014) was most likely due to assuming that the 
high— J CO lines originate in the dense gas traced 
by HCN, whereas in NGC 6240 they originate in 
shocks (Meijerink et al. 2013) and the HCN lines 
are cold, dense gas dominated, isolated by using 
our new H 13 CN observations. Our models suggest 
instead that [ 12 C]/[ 13 C] = 98^° and is “standard” 
for ULIRGs, and even consistent with the Milky 
Way average value of 68. The high [HCN]/[H 13 CN] 
and [HCO + ]/[H 13 CO + ] ratios (3002 qq) are due to 
isotope fractionation in cold, dense gas, with the 
[HCO + ]/[H 13 CO + ] ratio being particularly uncer¬ 
tain. 

• derived aco values for the diffuse, shocked and 
dense gas phases in NGC 6240, and combined these 
to produce a global aco = 1.5^ 4 , in between the 
ULIRG value of 0.8 and 4.2 appropriate for Milky 
Way GMCs. We also calculate ancN = 32f§, where 
«hcn is tracing the cold, dense gas. 

These results demonstrate the extreme value of 
isotopologue lines in LVG modelling, with their 
powerful ability to break degeneracies even when 
the isotopologue abundance ratios are unknown. 
Similarly, we have highlighted the dangers of fit¬ 
ting multiple species gas phase components by eye, 
as was quantified for multiple CO gas phases by 
Kamenetzky et al. (2014). 

Our model produces chemically feasible molecu¬ 
lar abundances and isotopologue abundance ra¬ 
tios, despite the significant simplifications. Nev¬ 
ertheless, the high [HCO + ]/[H 13 CO + ] ratio is un¬ 
expected, and warrants further investigation. A 
distinct possibility, evidenced by our slightly low 
continuum flux density, is that we are missing 
a significant fraction of the molecular line flux. 
This would present as elevated [HCN]/[H 13 CN] and 
[HCO + ]/[H 13 CO + ] abundance ratios. Importantly, 
even if this is the case the main conclusions of this 
work, that the elemental [ 12 C]/[ 13 C] ratio is not 
elevated in NGC 6240, are not affected; they are 
instead reinforced if we are missing line flux. 

In any case, further observations of 13 CO, H 13 CN 
and H 13 CO + are warranted, both in NGC 6240 
and other galaxies with CO, HCN and HCO + 
SLEDs. These isotopologue lines provide added 
value and, as we have demonstrated here, even 
a single line can simultaneously improve temper¬ 
ature constraints while constraining the isotopo¬ 
logue abundance ratio, without the need to assume 
a, potentially erroneous, ratio. 
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APPENDIX 
A. LVG MODELS 

We wrote an MCMC wrapper for the LVG code RADEX to allow a rapid and efficient exploration of high dimension 
parameter spaces and easy generation of posterior pdfs. The code allows us to fit multiple species simultaneously, as 
well as fitting multiple gas phases simultaneously. The code is written in Python, with the RADEX calculations in 
Fortran presenting the rate limiting step: > 99.8% of the runtime is RADEX IO and calculations. 

We adopt a Metropolis-Hastings (MH) sampler for its simplicity and robustness. The likelihoods from the LVG 
models are generally simple, contiguous shapes, and as such are well suited to the MH algorithm. We do not include 
covariance terms between the new parameters when sampling: each is sampled from a ID Gaussian distribution. We 
adopt a x 2 likelihood function: 


C oc exp (-X 2 / 2 ) ■ (Al) 

Constraints on the parameter ranges are implemented using uniform priors in log-space. The MCMC is conducted 
in log-space, with values converted to linear space for input to RADEX. This provides simple scale invariance and 
prevents biasing towards high values of the parameters. 

We use stock RADEX, and as such must be wary of non-convergences, bad fluxes and high optical depths. These are 
checked for when reading the RADEX results, and if identified the point is rejected and resampled. The alternative 
would be to treat this as a point of zero likelihood - test runs of these two showed no differences in the results. Since 
the MCMC naturally tends towards regions of greater likelihood the RADEX errors are not usually significant except 
for when using randomised initial positions. In the rare instance that the code becomes trapped in a region of errors 
new initial parameters are randomly generated. 


LVG Parameters 

Our aim was to extend the results of Papadopoulos et al. (2014) with our new data and an MCMC model. We 
therefore adopted a similar parameter space and K vir constraints: 0.5 < K vir < 20. See Table 3 for the three phase 
parameter ranges and Table 6 for the precursor model inputs. 

While we model the abundance and velocity gradient A' mo i and du/dr, RADEX requires as inputs the column density 
Nmol and line width Av. As only the ratio N mo i/Av is relevant 15 , we fix Av = lOkrns -1 and find N mo \ from: 


-Vnol 


3.08 x 10 18 7iH 2 X mo iAu 
du/d r 


(A2) 


B. CLOUD KINEMATIC SEPARATION 

Our MCMC LVG model fits line ratios, so while it is ideal for identifying the dominant conditions of the sampled 
gas phases it is poorly adapted for identifying the masses within each phase. A simple but inelegant solution would be 
to identify the CO line luminosity from each gas phase and then use a canonical aco factor to estimate the mass of 
the phase. This is however extremely susceptible to uncertainties in the aco factor. Instead, we derive an alternative, 
self-consistent, method, which is used to generate aco factors for each gas phase and for the galaxy as a whole. We 
apply a physically motivated upper limit on N^ 2 , thereby applying a limit on Au per cloud. 

A single gas cloud representing the galaxy is inadequate, since the column density of H 2 can be written as: 


n h 2 


3.08 x 10 18 ?r H2 Au 
clu/d r 


(Bl) 


For the high density gas phase this leads to column densities in excess of 10 26 cm -2 ! In reality, the dense gas is not in 
a single cloud but distributed across the galaxy in many, much smaller clouds, so that Au is determined by the large 
scale motion in the galaxy and does not correspond to the LVG Au, which is the line width of the putative emitting 
cloud. 


15 For a given product N mo i/ Av, changing Av only changes the 
quoted flux, proportional to the change in Au, which cancels in 


any ratios. The optical depths and brightness temperatures remain 
unchanged. 
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We therefore separate the emission into multiple identical clouds to obtain more realistic estimates of the mass of 
the gas phases, without assuming cloud line widths. We start with the general statement that for a cloud of uniform 
density: 


47t ( (jAv \ 

= T( v 2 duM ; J nM 


(B2) 


where ( = 3.08 x 10 18 is the parsec to cm conversion factor, Av is the FWHM line width of the cloud, du/d r is the 
velocity gradient of the cloud in kms _1 pc _1 , ?ih 2 is the hydrogen number density in cm~ 3 and n h 2 is the helium 
adjusted ISM mass per H 2 molecule. 


Separately, we make the general statement of beam dilution that for a Gaussian surface brightness with peak Tb and 
a Gaussian beam much larger than the size of the cloud: 


^bX 2 1-133 / Av A 2 

Ubeam \ 2(l l' / ( 1 f‘ J 


(B3) 


where x = ygy arcsec/pc is the angular scale for NGC6240 and flbeam is in arcsec 2 . In this case, Tb is a result of 
RADEX runs, which give optically thick clouds for the HCN and CO J = 1 — 0 lines. 


For a Gaussian line, or a line which can be described as a combination of Gaussians, the number of clouds is the ratio 
of the area (TAv) of the observed line to the diluted area per cloud: 


_ T 0 bsFWHM 

^clouds m \ 

T h ,d&v 


TobsflbeamFWHM 

1.133TbAvx 2 ( 23 ^ 73 ?) 2 


Therefore the total mass in the galaxy in this gas phase is: 


(B4) 


Af — A/cioud ^clouds 


(B5) 


4-7T / (Av 
3 \2dv/dr 


TobsObeamFWHM 

^H 2 MH 2 ~2 

1.133T b Avx 2 (^) 


1.85C 3 Mh 2 T obs O b earn ^H 2 
X 2 Tb dz> / d7' 


(B6) 


(B7) 


= 0 i2 3 ?obs^beam J/H2 _fWHM [Mq] (B8) 

TbX 2 clv/dr 

This in fact corresponds to the same equation for a mo i as that derived by Papadopoulos et al. (2012) by considering 
iso-velocity surfaces across a galaxy. Also, we find that the mass is independent of the chosen upper limit on 7 Vh 2 
which is reassuring given that GMCs are likely to present a range of column densities (Vazquez-Semadeni et al. 1997). 

C: PRECURSOR LVG MODELS 

The three phase model presented in Section 4.1 makes assumptions regarding the gas phases and distribution of 
molecular species between phases. These assumptions are based on a series of models of increasing complexity that 
are precursors to our three phase, six species model. These models do not add to the scientific results of the paper, 
but due to their role in justifying our assumptions we include them here. 

Single Phase 

We first ran four models with single species and their respective isotopologues: HCN+H 13 CN with the detected 
H 13 CN(1 - 0) line flux, HCN+H 13 CN with the upper limit on the H 13 CN(1 - 0) line flux, HCO++H 13 CO+ and CS 
only, with CS data shown in Table 5. We also run a combined model with HCN+HCO + +H 13 CN+H 13 CO + . The input 
parameter ranges for these models are shown in Table 6. When fixed, molecular abundances were chosen to match 
those of Papadopoulos et al. (2014), i.e., Ahcn = 2 x 10" 8 , X HCO + = 8 x 10 -9 and Acs = 1 x 10” 9 . 

The numerical results of the models are given in Table 7, where we present both the the means of the MCMC traces 
and the best fit values. We recover similar temperatures and densities to Papadopoulos et al. (2014), although the 
inclusion of our new 13 C isotopologue lines allows us to place tighter constraints on these parameters. 
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TABLE 5 

LVG OS Line Inputs 


Line 


^0 

[GHz] 

*$line 

[Jy kms -1 ] 

References 21 

CS(2 — 

1) 

97.98 

7.5 ± 1.5 

i 

CS(3 — 

2) 

145.0 

9± 2 

i 

CS(5 — 

4) 

239.1 

54 ±6 

2 

CS(7 — 

6) 

342.9 

7.4 ±1.0 

3 


a l = Papadopoulos et al. (2014), 2 = Wang et al. (2011), 3 = Scoville et al. (2015) 


TABLE 6 

MCMC MODEL INPUT PARAMETER RANGES. 


-Ik 

Model_ID [K] [cm -3 ] 


HCN & H 13 CN, measured. 

la 

0.5 - 3 

2-8 

HCN & H 13 CN, upper limit. 

lb 

0.5 - 3 

2-8 

HCO+ & H 13 CO+. 

2 

0.5 - 3 

2-8 

CS. 

5 

0.5 - 3 

2-8 

HCN, HCO+, H 13 CN & H 13 CO+ 

4 

0.5 - 3 

2-8 


dn/dr 

[kms -1 pc -1 ] 
0-3 
0-3 
0-3 
0-3 
0-3 


All values log 10 

[HC(N/Q+j] 
[Hl 3 C(N/0+)] 


A'hCN 


A'siO - Y HCO+ 


1 - 3 -7.70 

1 - 3 -7.70 

1-3 


- 8.10 


1 - 3 


-12-4 


-12-4 


Acs 


-9.00 


TABLE 7 
MCMC Results 


Posterior mean results and the upper and lower limits on the 68% credible interval. 


Model 

ID 

T k 

[K] 

nu 2 

[cm -3 ] 

dv/dr 

[kms -1 pc -1 

All values log 10 
[HC(N/0+)] 
[h 1 3c(n/o+)] 

AhCN WiO 

A'hCO+ 

HCN & H 13 CN, measured . 

HCN & H 13 CN, upper limit. 

HCO+ & H 13 CO+. 

CS. 

HCN, HCO+, H 13 CN & H 13 CO+ 

la 

lb 

2 

3 

4 

1 ' 52 L 00 

1 cnl.75 
i,OU 0.95 

1 oq 1.60 
1,oy 1.13 

2 71 30 

1 x 2.65 

1 f)4 1,13 
i - U4 0.85 

4 no5.60 
^• yo 4.40 

4 . 95 |;!f 

o o4.1 
°-°3.0 
c 076.62 
'5.42 

1 OR 2 - 28 
±,ZO 0.42 

1 oql.98 
1 - oy 0.42 

1 QQl-05 
1 - UU 0.03 
n «o0.93 
u -°°0.03 

1 qi2.34 

9 972.46 
^■^‘1.90 

9 403 .OO 
z '^ y 2.34 

9 743 .OO 
'^ 2.68 

9 ni 2 -82 
z,01 2.42 

0 09 — 7.52 

0 .o^_9 12 

00 
06 oi 
. 1 1 

: : : : ^ 
0 
ci 

1 


Joint maximum likelihood results. 

Model 

ID 

T k 

[K] 

n n 2 

[cm -3 ] 

dv/dr 

[kms -1 pc -1 

All values log 10 
[HC(N/0+)] 

AhCN 

Asio 

X HCO+ 

[H 13 C(N/0+)] 

HCN & H 13 CN, measured . 

la 

0.86 

6.3 

2.7 

3.0 




HCN & H 13 CN, upper limit . 

lb 

0.88 

6.3 

2.7 

3.0 




HCO+ & H 13 CO+. 

2 

1.1 

4.8 

1.1 

3.0 




CS. 

3 

3.0 

3.5 

0.1 





HCN, HCO+, H 13 CN & H 13 CO+ 

4 

0.96 

6.3 

1.7 

2.6 

-9.0 


-9.8 


The Effect of T bg 

For some very active galaxies with luminous dust emission, background radiation from hot dust behind the gas cloud 
may significantly affect the line ratios (Papadopoulos et al. 2010). For dense molecular gas tracers, neglecting this 
may lead to extremely unusual solutions from RADEX fits to the line ratios. 

We explored this possibility in NGC6240 by testing whether increasing T bg can improve the fits of our HCN + 
HCO + single phase model. The model was rerun with T bg = 10, 20 and 50 K. 

We find that increasing T bg to 10 K marginally improves the HCO + SLED fit at the expense of a sightly worse HCN 
SLED fit. The solutions are constrained to give positive fluxes which indirectly constrains the gas kinetic temperature 
in these models to be > T bg . The results push against this limit, so that Tk/Tiust ~ Tk/Tbg — 1, which is unlikely to 
be reasonable in the extensively shocked environment of NGC6240 (Meijerink et al. 2013; Papadopoulos et al. 2014), 
unless the HCN and HCO + are located not in shocks but in cold, dense cores. As T bg is increased, the fits become 
progressively worse, with 50K not remotely fitting the SLED and producing \ 2 > 1300 (cf. 4 — 6 for T bg = 3 — 10K). 
The SLED fits for these models are shown in Figure 13. Note that the change in shape is largely due to a suppression 
of the low—J line emission in excess of the background due to the increased background. 

The 50K case is comparable to the the 52K black body dust temperature described by Tacconi et al. (1999). 
However, we found no consistent solutions over multiple runs, and pathological fits to the SLEDs for this T bg . This 
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Fig. 13.— SLED fits for the HCN+HCO"*” single phase model, with Tb g increasing as 3, 10, 20 and 50 K. The best fits are from the 3K 
background, and the two observationally motivated background temperatures (20 and 50 K) both produce much worse fits. 

is consistent with the inconsistencies between the CO and 1.3 mm continuum estimates of gas mass commented on 
by Tacconi et al. (1999): they find that the dust co-located with the dense gas was most likely significantly cooler 
(perhaps 20-25 K), and furthermore that the shock heated gas is thermally decoupled from the cooler dust. 

The inability of our results to find good a solution for a 50 K background is unsurprising, but the poor solution at 
20 K is slightly more interesting. One possible, and we would argue likely, explanation could simply be that the HCN 
and HCO + emission are tracing the relatively quiescent regions of shielded, dense gas, i.e., they are performing as 
classical dense gas tracers. This is consistent with the marginal improvement (worsening) of the SLED fit of HCO + 
(HCN) when Tb g is increased from 3K to 10 K, as HCO + is probably more extended and less shielded than HCN, so 
subject to a slightly stronger background radiation field. It is also consistent with the very high HCN/CN line ratio 
of 2 found by Aalto et al. (2002), which is evidence of large reservoirs of UV-shielded gas (Fuente et al. 1995; Greaves 
& Church 1996). 

We conclude that there is insufficient evidence to support using a super-CMB background and proceed with the 3 K 
CMB. 


Two Phase Modelling - HCN+HCO + 

Here, we extend our model to include a second, non-interacting gas phase with both phases containing HCN and 
HCO + . This is strongly physically motivated as single dish observations across the galaxy are sampling multiple 
non-interacting gas phases confused within the beam. 

Each of the two phases has their own kinetic temperature, density, and velocity gradient while they share common 
HCN and HCO + abundances and a single [HC(N/0 + )]/[H 13 C(N/0 + )] abundance ratio 16 . The fluxes from the two 
phases are added, weighted by a free parameter, /, such that: 

S = S 1 f + S 2 (l-f), (Cl) 

where Si and S 2 are the RADEX fluxes for phase 1 and 2 respectively, in a simplified version of the method used in 
(Zhang et al. 2014). Either of Si and S 2 may be negative, but the combined fluxes S are constrained to be positive, 
as per observations. Furthermore, there is no line transfer between the two phases: they emulate two phases confused 
in the imaging beam, not two phases along the same line of sight (unless they are decoupled by a sufficiently large 
velocity offset). 


The results of this two phase model are shown in Figure 14. It appears from Figure 14 that we have an underdeter¬ 
mined model unable to either constrain both phases or clearly exclude one. Even though we are unable to constrain 
the kinetic temperature, density or velocity gradient, the abundance ratio is well determined. 

The [HC(N/0 + )]/[H 13 C(N/0 + )] ratio in this model is elevated. It is almost definitely unsuitable to have a single 
[HC(N/0 + )]/[H 13 C(N/0 + )] ratio for both phases and both HCN and HCO + ; unfortunately, as is shown by the 


16 While this is undesirable it is necessary to minimise the num¬ 
ber of free parameters we attempt to constrain with our few line 


ratios. A model was attempted with free abundances in each phase, 
but this was unable to converge on any solutions. 
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Fig. 14.— The modified step plot for the two phase HCN+HCO”*” model. Contours are at the 68%, 95% and 99% credible intervals. 
Red/solid and blue/dashed contours are phases 1 and 2, with Tbg — 3 and 3K respectively. [ 12 C]/[ 13 C], Xhcn and ^HCO+ are kept the 
same in the two phases, so have only one set of black contours. The two phases are completely degenerate - unsurprising given the good 
fits of the single phase models. 

degeneracy of the two phases, we do not have sufficient observations nor are observations of sufficient precision to 
constrain two phases without making some simplifications. Further observations of the H 13 CN and H 13 CO + SLEDs, 
in particular the J = 3 — 2 line, would go a long way towards solving this issue. 

Two Phase Modelling - HCN+HCO + +CO 

The previous two phase model failed due to the paucity of available molecular lines and the large uncertainties on 
those available. In lieu of additional 13 C isotopologue line observations, which would place much tighter constraints 
on the models, we explore the addition of CO and 13 CO to the two phase models. We take the line fluxes directly from 
Papadopoulos et al. (2014), and do not make any corrections for source size: the / scaling parameter allows the model 
to account for different beam filling factors due to the different sizes of the HCN and CO emitting regions, while the 
two phases allow for non-cospatial emission. 
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TABLE 8 

HCN+HCO++CO Two Phase Model Solutions. All values log 10 . 


Parameter 

T k ,i. 
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-5-3=e:o 
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-O-Hais 

-5.2 
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With the additional lines we are able to extend the model complexity so that the free parameters include Xqo , 
Xhcni ^hco+ and [Xi2c]/[Xi3c] in each gas phase. Once again, we have insufficient lines to simultaneously fit free 
X mo i and [Xl2c]/[Xi3c] ratios for every species in both phases, and we use a single [Xi 2 c]/[Xisc] ratio in each phase 
for all species. Even with these restrictions the code can take many more steps to converge, and the traces were 
carefully inspected to exclude unstable/unconverged runs. Of the runs that converged, all converged onto the same 
locus, which is presented here. The results are tabulated in Table 8 and shown in Figures 15 (parameter pdfs) and 16 
(SLEDs). 

This model provides very tight constraints on the two gas phases, although the tightness of the constraints is no 
guarantee of accuracy. The model spontaneously separates into a hot, diffuse, CO dominated phase and a cold, dense, 
HCN dominated phase. This result is robust to randomised initial positions over multiple MCMC runs and is similar 
to the findings of Kamenetzky et al. (2014) who used only CO. T k and nn 2 are the best constrained parameters, with 
less than ±0.1 dex spread in the 68% credible intervals. In the diffuse phase Xhcn falls rapidly, and Figure 16 confirms 
that HCN has no contribution from the diffuse gas phase, so that Xhcn in this phase is an upper limit only. In both 
phases Xco is slightly lower than the canonical value of 2 x 10~ 4 presenting instead 4 x 10~ 5 and 6 x 10 -6 . However, 
the spread in Xco is very large, especially in the hot phase. 


The [Xi2 G ]/[Xi3 G ] ratio: evidence for fractionation? 

Even with the addition of the CO and 13 CO lines we are unable to constrain two gas phases with free CO, HCN 
and HCO + ; 13 CO, H 13 CN and H 13 CO + abundances in each phase. In the quest to identify whether the [ 12 C]/[ 13 C] 
ratio found by Papadopoulos et al. (2014) is real, or due to biased SLED fitting, we ran a model with Xco, X H cn 
and X HC o+ fixed to the best fit parameters of the previous run (see Table 8 for values). We then introduce as free 
parameters [CO]/[ 13 CO]i, [CO]/[ 13 CO] 2 , [HCN]/[H 13 CN] and [HCO+]/[H 13 CO+]. I.e., while HCN and HCO+ are 
present in both phases, each has a single isotopologue ratio used for both phases. This is motivated by us only having 
the J = 1 — 0 line for the H 13 CN and H 13 CO + , which have negligible contributions in the hot phase - in effect we are 
indirectly excluding H 13 CN and H 13 CO + from the hot phase. 


The results of the [ 12 C]/[ 13 C] investigation model are shown in Figure 17 and Table 9. The model recovers very 
similar T k , nn 2 and du/dr as the previous model. As we have fixed the main isotope abundances this is not surprising, 
but it is an important confirmation that the isotopologue abundance ratios we derive here are compatible with the 
less restricted two phase model; by fixing some parameters and freeing others we have not inadvertently moved into a 
very different locus in parameter space. 

The [HCN]/[H 13 CN] and [HCO+]/[H 13 CO+] ratios are both very high, while the [CO]/[ 13 CO] ratio is ~ 100 - 200 
in both phases. However, the abundance ratios are correlated, so that the spread in ([CO]/[ 13 CO] 2 )/([HCN]/[H 13 CN]) 
is extremely small, strongly peaking at 2 — 3, consistent with Roueff et al. (2015). These results should be interpreted 
with caution: the fixing of the molecular abundances could be biasing the results. However, it is expected to bias 
them towards the values found in the previous model, so [CO]/[ 13 CO] 2 lying much lower, while [HCN]/[H 13 CN] and 
[HCO + ]/[H 13 CO + ] both lie higher, suggests that this is a real effect. This is evidence for a “standard” [ 12 C]/[ 13 C] 
abundance ratio (for (U)LIRGs) of ~ 100 in NGC6240. The high ratio found by Papadopoulos et al. (2014) is 
explained by their attribution of the high—J CO lines to the dense gas phase, whereas our simultaneous fitting and 
the addition of the H 13 CN and H 13 CO + lines/limits leads to the model spontaneously fitting these high—J lines to a 
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Fig. 15.-— The modified step plot for the two phase HCN+HCO~*~+CO model. Contours are at the 68%, 95% and 99% credible intervals. 
Red/solid and blue/dashed contours are phases 1 and 2, with Tb g — 3 and 3K respectively. Contours of constant pressure have been 
included in the — nn 2 plot at intervals of one dex. There is a clear, robust and spontaneous separation into a hot, diffuse and cold, 
dense gas phase. 

hot, diffuse phase, as did the two phase CO model of Kamenetzky et al. (2014). Our high observed [HCN]/[H 13 CN] 
and [HCO + ]/[H 13 CO + ] ratios are explained by chemical effects (ICE) in the very cold, dense gas. 
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Fig. 16.— SLED fits for the two phase HCN+HCO^+CO model, showing the relative contribution of each gas phase. The same 
normalisation has been used for the CO and 13 CO SLEDs. 
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Fig. 17. — ID marginalised posterior pdfs for the gas phase parameters of the [Xi 2 q\/[Xizq\ investigatory model. The gas parameters 
are very similar to those of the HCN+HCO”*~+CO 2 phase model. The [CO]/[ 13 CO] abundance ratios are similar in both phases, and 
are in the “standard” region for (U)LIRGs, suggesting that the high [HCN]/[H 13 CN] and [HCO+J/fH-^CO -1- ] ratios are due to isotope 
fractionation in the cold gas phase. 


TABLE 9 

[ 12 C]/[ 13 C] Study Two Phase Model Solutions. All values log 10 . 
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